{ "cells": [ { "cell_type": "markdown", "id": "endangered-willow", "metadata": {}, "source": [ "# Weighted Nonlinear Fit\n", "*March 18, 2021*\n", "\n", "In this Jupyter notebook we will fit a nonlinear function to a set of experimental data. The fit will be weighted by the measurement uncertainties.\n", "\n", "In this example, we take data from a series $LRC$ circuit. For the $y$-axis we will have the ratio of the magntidue of the voltage across the resistance to that supplied by the function generator driving the circuit ($V_\\mathrm{ratio}=V_R/V_0$). On the x-axis we will have angular frequency $\\omega$.\n", "\n", "The process will be very similar to the process used to fit a linear function to a dataset. The main difference will be in the definition of the model function that we are fitting to.\n", "\n", "- First, import the *NumPy*, *Matplotlib*, & *SciPy* modules." ] }, { "cell_type": "code", "execution_count": 16, "id": "incident-korean", "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt\n", "from scipy.optimize import curve_fit" ] }, { "cell_type": "markdown", "id": "entitled-dominant", "metadata": {}, "source": [ "- Next, enter the data (frequency and $V_\\mathrm{ratio}\\pm\\Delta v_\\mathrm{ratio}$) as arrays." ] }, { "cell_type": "code", "execution_count": 17, "id": "chronic-pricing", "metadata": {}, "outputs": [], "source": [ "Vratio= np.array([0.198e-1, 0.67e-1, .117, .185, .331, .450, .573, .689, .718,\\\n", " .714, .704, .670, .631, .549, .412, .318, .249, .186, .138])\n", "errVratio= np.array([0.1e-2, 0.2e-2, 0.3e-2, 0.3e-2, 0.5e-2, 0.6e-2, 0.7e-2,\\\n", " 0.8e-2, 0.8e-2, 0.8e-2, 0.8e-2, 0.8e-2, 0.8e-2, 0.7e-2,\\\n", " 0.6e-2, 0.4e-2, 0.5e-2, 0.3e-2, 0.3e-2])\n", "f = np.array([3000, 10000, 15000, 20000, 25000, 28000, 30000, 32000, 33000,\\\n", " 34000, 35000, 36000, 37000, 39000, 43000, 47000, 52000, 60000,\\\n", " 70000])" ] }, { "cell_type": "markdown", "id": "veterinary-participation", "metadata": {}, "source": [ "- Calculate the angular frequency $\\omega$." ] }, { "cell_type": "code", "execution_count": 18, "id": "talented-trailer", "metadata": {}, "outputs": [], "source": [ "omega = (2*np.pi*f)" ] }, { "cell_type": "markdown", "id": "imposed-basketball", "metadata": {}, "source": [ "- Plot the data using *errorbar(x, y, e)*." ] }, { "cell_type": "code", "execution_count": 19, "id": "municipal-chinese", "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEGCAYAAABo25JHAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAhZ0lEQVR4nO3de3xcdZ3/8debSbGCWy5SFdpAC8vFwk/SEltRl4XKpYXUyIIPuVSrrttfFVZdcRV+Ig3t+sPL6k9ccWth2dUtFy9cLJG2SBVdL902kICltFi5NUVpuQgIujTh8/vjnNBJOplMmjmZTOb9fDzmkTnfc+bMZ75N85lzvud8vooIzMysdu1R6QDMzKyynAjMzGqcE4GZWY1zIjAzq3FOBGZmNa6u0gEM1gEHHBCTJk2qdBhmZlXl7rvvfjIixhdaV3WJYNKkSbS1tVU6DDOzqiLp0f7W+dSQmVmNcyIwM6txTgRmZjXOicDMrMY5EZiZ1TgnAjOzGudEYGZW46ruPgKzLHR3w4oV0N4OU6fC7NmQy1U6KrPh4SMCq1ktLS1IQsoxdtxK5l20lRXP38u8i7YydtxKpBySaGlpqXSoZplStU1M09jYGL6z2MqptRUuWPgsc9csJzcm6N4h/nXKTJ7ZfA3QBbQDK4CXAVi4cKGTg1UdSXdHRGOhdT41ZDWvvR3qT91CbkzypUh7QG7M3rz2iI9wZPPv6Fz9SaZM2JvWW+p8ushGJScCq3lTp8I1C+vpXtRObkzw4G0TIMT/Xr8yPUK4l2UzmmltHUcu53EEG30yHSOQNEvSJkmbJV1cYP0/SupIH+sldUvaP8uYzPqaPRumTNibbx5zMqs/1cAdH3krR7xz5xFCbkxw0MyHOfu8TR5HsFEpszECSTngQeAUoBNYB5wbERv62X4O8A8RMbPYfj1GYOXS0tLC5Zdfni7tAcwGGoAx7H/Eh1iw/sevjBlcfdRZ7PmaHcxr2zmOsGxGM1ctGkdTU8U+glnJio0RZHlEMB3YHBEPRcRLwI1Ac5HtzwVuyDAes15aWlqIiPTRzcKFjcDngEU81/lrlhw9kzs/dSxLjp7JU1u2cvBpj/Q6Sqg/7TE6Oir5CczKI8tEMAHYkrfcmbbtQtJewCzgpgzjMStqZ2Lo5s/PzeJbX5nA6eOO5VtfmcDN3zucrasn0b1DAHTvEFtWHUxDQ2VjNiuHLBOBCrT1dx5qDvCLiHi64I6k+ZLaJLVt3769bAFa7dh5z4DS8/pNSJ9Nf+ZeWddzrj+Xg6YmuPTS5GdTU+9xhGUzmpkycS9mz67s5zIrhyzHCI4HWiLitHT5EoCIuKLAtrcA34uI6wfar8cIbCi6u6HpzC7WPbiNI9/5JJ2rDyl6aWj/4wgd+N4CqybFxgiyTAR1JIPF7wC2kgwWnxcR9/fZbh/gYaA+Il4YaL9OBDYUhW4e86Cv1YKKDBZHRBdwIbAKeAD4bkTcL2mBpAV5m54J3FFKEjAbqr43j3nQ1yzj+wgi4vaIOCIiDouIz6VtSyJiSd42/xER52QZh1mPqVNhyx31HvQ1y+Oic1YTegaL58zJ8fjGX/W6NPTxTb9kzpzB3xjW3Z2calq8OPnZ3Z1d/GZZctE5qzk9Jac7OqChYXClInYOHu9B3V4/ZNzE/8WRzU+y6QcH8Fznr+l68QzgZQ8e24hTkcHirDgR2EjgQWerNpW6s9hs1PKgs40mTgRmu8GDzjaaOBGYDUIWg85mleYxArPdNJRBZ7Ph5hnKzDLQU4/Ig8NW7XxqyMysxjkRmJnVOCcCM7Ma5zECszLqGUD2BPdWTXxEYDZEOye9yTF23EpPcG9Vx5ePmpWJy07YSOYSE2bDwGUnrFo5EZiVictOWLVyIjAbIpedsGrnMQKzMnLZCRupXGLCbJi47IRVo0xPDUmaJWmTpM2SLu5nmxMldUi6X9JPs4zHzMx2ldkRgaQccBVwCtAJrJO0PCI25G2zL/ANYFZEPCbpdVnFY2ZmhWV5RDAd2BwRD0XES8CNQHOfbc4Dbo6IxwAiYluG8ZiZWQFZJoIJwJa85c60Ld8RwH6S7pJ0t6T3FdqRpPmS2iS1bd++PaNwzcxqU5aJQAXa+l6iVAccB5wBnAZ8VtIRu7woYmlENEZE4/jx48sfqZlZDcvyqqFOoD5veSLweIFtnoyIF4AXJP0MOBZ4MMO4zMwsT5ZHBOuAwyVNlrQncA6wvM82PwD+SlKdpL2AGcADGcZkZmZ9ZHZEEBFdki4EVgE54NqIuF/SgnT9koh4QNJK4D7gZeCaiFifVUxmZrYr31lso4rnAzArzNVHbVTzfABmQ+MjAhs1PB+AWf98RGA1wfMBmO0eJwIbNTwfgNnucSKwquf5AMyGxmMENqp4PgCzwjwfgdWMapgPwJe42kjjU0Nmw8CXuNpI5lNDZsPIl7hapfjyUbMRwpe42kjkRGA2jHyJq41ETgRmw8CXuNpI5jECs2HmS1ytEnz5qNkIUg2XuFpt8akhM7Ma50RgZlbjnAjMzGqcE4GZWY1zIjAzq3GZJgJJsyRtkrRZ0sUF1p8o6VlJHenjsizjMTOzXWV2+aikHHAVcArQCayTtDwiNvTZ9L8iwhfSmZlVSJZHBNOBzRHxUES8BNwINGf4fmZmthuyTAQTgC15y51pW1/HS7pX0gpJRxfakaT5ktoktW3fvj2LWM3MalaWiUAF2vrWs7gHOCQijgX+Bbi10I4iYmlENEZE4/jx48sbpZlZjcsyEXQC9XnLE4HH8zeIiOci4o/p89uBMZIOyDAmMzPrI8tEsA44XNJkSXsC5wDL8zeQ9AZJSp9PT+N5KsOYzMysj8yuGoqILkkXAquAHHBtRNwvaUG6fglwNvBhSV3An4BzotrKoZqZVbmSylBLeidwQrr404i4LdOoinAZajOzwRvSVJWSrgA+BmxIHx9N28zMbBQo5dTQGUBDRLwMIOlbQDtwSZaBmZnZ8Ch1sHjfvOf7ZBCHmZlVSClHBFcA7ZJ+QnJvwAn4aMDMbNQYMBFExA2S7gLeTJIIPh0Rv886MDMzGx79JgJJR0XERknT0qbO9OdBkg6KiHuyD89sp55J39vbYepUT/puVi7Fjgg+AcwHvlxgXQAzM4nIrIDubmg6s4sNW1+g/tQtXLOwnilL96b1ljonA7Mh6newOCLmp09nR8RJ+Q/g9OEJz2pdS0sLkqira2LtpieYu2Y5M6+4h7lrlrN24xPU1TUhiZaWlkqHala1Srlq6JcltpmVXUtLCxHBokWtvPFvniE3JrkBMjcmeONZz7B4cSsR4URgNgT9JoK0DtBxwKslTZU0LX2cCOw1XAGaQTImsOWOerp3JEVtu3eILasOpqGhsnGZjQbFjghOA/6ZpGroV0jGCr5MMnbwf7IPzWznqaE5c3I8vvFXLDl6Jnd+6liWHD2Txzf9kjlzcj41ZDZEA9YaknRWRNw0TPEMyLWGalfPVUMdHdDQ4KuGzAajWK2hUu4juEnSGcDRwNi89kXlC9FsYLkcNDUlDzMrn1KKzi0B3gP8PckNZe8GDsk4LjMzGyalXDX01oh4H/BMRFwOHE/vmcfMzKyKlZII/pT+fFHSQcAOYHJ2IZmZ2XAqpehcq6R9gS+RTDYfwDVZBmVmA3PJDSuXAY8IImJxRPwhvXLoEOAo4POZR2Zmu+i5nFbKMXbcSuZdtJUVz9/LvIu2MnbcSiRfTmuDV/TyUUkTgAOB+yLiJUmvAz4OvD8iDhqeEHvz5aNm0NoKFyx8lrlrlpMbE3TvEMtmNHPVonG+qsoK2q2pKiV9HOgA/gVYI2ke8ADwauC4Et94lqRNkjZLurjIdm+W1C3p7FL2a1br2tuh/tQtvUpu1J/2GB0dlY3LqlOxU0PzgSMj4njgXcDVwBkR8Q8R8buBdiwpB1wFzAamAOdKmtLPdl8AVg0+fLPa5JIbVk7FEsGfI+JpgIh4DHgwItYMYt/Tgc0R8VBEvATcCDQX2O7vgZuAbYPYt1lNcskNy0K/YwSStpH88e5xTv5yRHy06I6T0zyzIuJD6fJ7gRkRcWHeNhOA60nmNvg3oDUivl9gX/NJjlA4+OCDj3v00UdL+nBmo5lLbthg7G6JiX/ss3z3YN+3QFvfrPNVkqkvu6VCm6cvilgKLIVksHiQcZiNSi65YeXSbyKIiG8Ncd+d9L4DeSLweJ9tGoEb0yRwAHC6pK6IuHWI721mZiUq5Yay3bUOOFzSZGAryaml8/I3iIhX7lCW9B8kp4ZuzTAmMzPrI7NEEBFdki4kuRooB1wbEfdLWpCuX5LVe5uZWemyPCIgIm4Hbu/TVjABRMT7s4zFzMwKK6UM9RGSVktany6/SdKl2YdmZmbDoZTqo1cDl5BUHSUi7iM5329mZqNAKYlgr4hY26etK4tgzMxs+JWSCJ6UdBjpPQDpjWIDlpgwM7PqUMpg8QUkN3MdJWkr8DAwN9OozMxs2JQyef1DwMmS9gb2iIjnsw/LzMyGy4CJQNIn+iwDPAvcHREd2YRlZmbDpZQxgkZgATAhfcwHTgSulvSp7EIzM7PhUMoYwWuBaRHxRwBJC4HvAyeQFKL7YnbhmZlZ1ko5IjgYeClveQdwSET8CfifTKIyM7NhU8oRwfUkU1X+IF2eA9yQDh5vyCwyMzMbFqVcNbRY0grgbSRzDCyIiJ7Z48/PMjgbnXomVGlvT6Zc9IQqZpVVyqkh0j/8NwA3A9skHZxpVDbq9EyxKOUYO24l8y7ayorn72XeRVsZO24lkqdYNKuUfqeqfGUD6Z3Al4GDSOYVPhjYGBFHZx/erhobG6OtrW3gDW1Eam2FCxY+y9w1y8mNCbp3iGUzmrlq0TjPtGWWoWJTVZZyRLAYeAvJ5PWTgZOBX5QxPqsh7e1Qf+oWcmOSLyC5MUH9aY/R0VHZuGpZd3eSoBcvTn52d1c6IhtupSSCHRHxFLCHpD0i4idAQ7Zh2WjTc2rossuaeODm/ejekcxR3b1DPHDTfnz2s00+NVQB3d3QdGYXFyx8llUvrueChc/SdGaXk0GNKSUR/EHSa4CfAddJuhJXH7VBamlpISLo6mpl+pGvZ9mMZn58yTSWzWhm+lGvp6urlYhwIhgmPYm5rq6JtZueYO6a5cy84h7mrlnO2o1PUFfnxFxLShkj2Bv4E0nSOB/YB1gWEU9nH96uPEZQ/XquGurogIYGXzVUSYsXw6oX1zPzinteafvxJdOYtfcxXOrpp0aVoY4RXBYRL0dEV0R8KyK+Bny6vCFaLcnloKkJLr00+ekkMPx8qs7ylXJEcE9ETOvTdl9EvGnAnUuzgCtJJq+/JiI+32d9M8lg9Mskp5s+HhE/L7ZPHxGYlU/PGMGGzhepP+0xtqw6mCkT96L1ljon6FGm2BFBvzeUSfow8BHgUEn35a36C0q4akhSDrgKOAXoBNZJWh4R+XcjrwaWR0RIehPwXeCogfZtZuWRy0HrLXWsWDGOjo5jaFjkU3W1qNidxdcDK4ArgIvz2p8vcXxgOrA5nc8ASTcCzeSVpegpZJfam3QWNDMbPj2n6nwfR+0qNkaQA54jmaHs+bwHkvYvYd8TgC15y51pWy+SzpS0Efgh8MFCO5I0X1KbpLbt27eX8NZmZlaqYkcEd7PzG7r6rAvg0AH23fc1Pa/r3RBxC3CLpBNIxgtOLrDNUpLpMmlsbPRRg5lZGfWbCNK7iIeiE6jPW54IPF7k/X4m6TBJB0TEk0N8bzMzK1EpZah76g2dkC7eFRGtJbxsHXC4pMnAVuAc4Lw++/1L4LfpYPE0YE/gqVKDNzOzoStlzuLPA28GrkubPibpbRFxSbHXRUSXpAuBVSTjDddGxP2SFqTrlwBnAe+TtIPkprX3xEDXs5qZWVmVch/BfUBDRLycLueA9lLuI8iC7yMwMxu8od5ZDLBv3vN9hhyRmZmNGKWMEVwBtEv6CcmVQCcARU8LmZlZ9Sh2Z/HXgesj4gZJd5GMEwj4dET8fpjiMzOzjBU7IvgN8GVJBwLfAW6IiI5hicrMzIZNv2MEEXFlRBwP/DXwNPDvkh6QdJmkI4YtQjMzy9SAg8UR8WhEfCEippLcB3Am8EDmkZmZ2bAo5T6CMcAskhvC3gH8FLg847jMbBTomYSovR2mTnVl05Gq3yMCSadIupakVMR84HbgsIh4T0TcOkzxmVmV6Zn0RsoxdtxK5l20lRXP38u8i7YydtxKpJwnvRlh+r2hLL1c9HrgpkpNS1mIbygzqw6trXDBwmeZu2Y5uTFB9w6xbEYzVy0a55LXFbBbN5RFxEkRcfVISgJmVj3a26H+1C3kxiRfNnNjgvrTHqOjo7Jx2a5KvbPYzKwkng+5+gxYa2ik8akhs+rg+ZBHlt2as9jMbCg8H3L1cCIws8x4PuTq4DECM7Ma50RgZlbjnAjMzGqcE4GZWY1zIjAzq3GZJgJJsyRtkrRZ0sUF1p8v6b708UtJx2YZj5mZ7SqzRJBOcn8VMBuYApwraUqfzR4G/joi3gQsBpZmFY+ZmRWW5RHBdGBzRDwUES8BNwLN+RtExC8j4pl0cQ0wMcN4zMysgCwTwQRgS95yZ9rWn78FVhRaIWm+pDZJbdu3by9jiGZmluWdxSrQVrCwkaSTSBLB2wutj4ilpKeNGhsbq6s4UhXzpCI22vh3urAsE0EnUJ+3PBF4vO9Gkt4EXAPMjoinMozHBuGVgmFbX6D+1C1cs7CeKUv3dsEwq1r+ne5flqeG1gGHS5osaU+SqS6X528g6WDgZuC9EfFghrHYIK1YARu2vsDcNcuZecU9zF2znA2dL7Ki4Mk7s5Grpyx2XV0Tazc90et3eu3GJ6irc1nszBJBRHQBFwKrSCa7/25E3C9pgaQF6WaXAa8FviGpQ5LrS1dYz3+aOXM+y4R3PNJrUpEJJz/MnDmX1vx/GqsuLS0tRASLFrXyxr95ptfv9BvPeobFi1uJiJr+nc70PoKIuD0ijoiIwyLic2nbkohYkj7/UETsFxEN6aNgrWyrhHY2/eCAXpOKbLr1AKCjolGZDZYnyhmYJ6axgjypiI02tf47XWxiGicC61fPFRYdHdDQ4CssrPrV8u+0E4GZWY0rlghcdM7MrMY5EZiZ1TgnAjOzGudEYGZW45wIzMxqnBOBmVmNy7LonJmZlUHWVVOdCMzMRrDhqJrqU0NmZiPQcFZNdSIwMxvRpnJk85O9qqYe+a4ngYayvYMTgZnZCNRTPvu22xazdfWkXlVTt945mdtu+6eylc/2GEGV89R7ZqPb7NkwZeneLJvR3Ktq6uzZ5XsPF52rYn0HkbbcUc+UCZ56z2y0KUfV1GJF53xEUMXyp5PMjQm6F7WzbEYzK1aMo6mp0tGZWbnkctDURGb/rz1GUIU8naSZlZMTQVXzdJJmNnSZjhFImgVcCeSAayLi833WHwX8OzAN+ExE/PNA+/QYwU61PvWemZWuImMEknLAVcApQCewTtLyiNiQt9nTwEeBd2UVx2iWy0HrLXWsWDGOjo5jaFjkq4bMbPCyHCyeDmyOiIcAJN0INAOvJIKI2AZsk3RGhnGMalkPIpnZ6JflGMEEYEvecmfaNmiS5ktqk9S2ffv2sgRnZmaJLBOBCrTt1oBERCyNiMaIaBw/fvwQwzIzs3xZJoJOoD5veSLweIbvZ2ZmuyHLRLAOOFzSZEl7AucAyzN8PzMz2w2ZDRZHRJekC4FVJJePXhsR90takK5fIukNQBswDnhZ0seBKRHxXFZxmZlZb5mWmIiI24Hb+7QtyXv+e5JTRmZmViGuNTQMXCHUzEYyJ4KMDcc0c2ZmQ+FaQxnLrxDaM83chs4XWbGi0pGZmSWcCDLiCqFmVi2cCDLnCqFmNrJ5hrKMuUKomY0EnqGsglwh1MxGOieCYeAKoWY2knmMwMysxjkRmJnVOCcCM7Ma50RgZlbjnAjMzGqcE4GZWY2rictHXf3TzKx/oz4RuPqnmVlxo/bUUE/Rt7q6JtZueqJX9c+1G5+grq7JRd/MzBjFiWCnqRzZ/GSv6p9HvutJoKGiUZmZjRSjNhG0tLQQEdx222K2rp7Uq/rn1jsnc9tt/0RE+IjAzGpepolA0ixJmyRtlnRxgfWS9LV0/X2SppU7htmzYcqEvVk2o5kfXzKNZTOamTJxL2bPLvc7mZlVp8wGiyXlgKuAU4BOYJ2k5RGxIW+z2cDh6WMG8K/pz7Jx9U8zs+KyvGpoOrA5Ih4CkHQj0AzkJ4Jm4NuRTIqwRtK+kg6MiN+VMxBX/zQz61+Wp4YmAFvyljvTtsFug6T5ktoktW3fvr3sgZqZ1bIsE4EKtPWdDq2UbYiIpRHRGBGN48ePL0twZmaWyDIRdAL1ecsTgcd3YxszM8tQlolgHXC4pMmS9gTOAZb32WY58L706qG3AM+We3zAzMyKy2ywOCK6JF0IrAJywLURcb+kBen6JcDtwOnAZuBF4ANZxWNmZoUpuWCnekjaDjwKHAA8WeFwRjL3T3Hun4G5j4qrtv45JCIKDrJWXSLoIaktIhorHcdI5f4pzv0zMPdRcaOpf0ZtiQkzMyuNE4GZWY2r5kSwtNIBjHDun+LcPwNzHxU3avqnascIzMysPKr5iMDMzMrAicDMrMZVXSIYaI6DaifpWknbJK3Pa9tf0o8k/Sb9uV/eukvSvtgk6bS89uMk/Tpd9zVJSttfJek7aft/S5qU95p56Xv8RtK8YfrIgyKpXtJPJD0g6X5JH0vb3UeApLGS1kq6N+2fy9N2908eSTlJ7ZJa0+Xa7p+IqJoHyR3KvwUOBfYE7gWmVDquMn/GE4BpwPq8ti8CF6fPLwa+kD6fkvbBq4DJad/k0nVrgeNJCvutAGan7R8BlqTPzwG+kz7fH3go/blf+ny/SvdHgf45EJiWPv8L4MG0H9xHSYwCXpM+HwP8N/AW988u/fQJ4HqgNV2u6f6peACD/Mc7HliVt3wJcEml48rgc06idyLYBByYPj8Q2FTo85OU8zg+3WZjXvu5wDfzt0mf15HcGan8bdJ13wTOrXRflNBXPyCZ/Mh9tGvf7AXcQzLZk/tnZ1wTgdXATHYmgprun2o7NVTS/AWj0OsjLcaX/nxd2t5ff0xIn/dt7/WaiOgCngVeW2RfI1Z6yD2V5Fuv+yiVnvboALYBP4oI909vXwU+Bbyc11bT/VNtiaCk+QtqSH/9Uayfduc1I46k1wA3AR+PiOeKbVqgbVT3UUR0R0QDyTff6ZKOKbJ5TfWPpCZgW0TcXepLCrSNuv6ptkRQq/MXPCHpQID057a0vb/+6Eyf923v9RpJdcA+wNNF9jXiSBpDkgSui4ib02b3UR8R8QfgLmAW7p8ebwPeKekR4EZgpqRl1Hr/VPrc1CDP7dWRDLBMZudg8dGVjiuDzzmJ3mMEX6L3QNYX0+dH03sg6yF2DmStIxkk7BnIOj1tv4DeA1nfTZ/vDzxMMoi1X/p8/0r3RYG+EfBt4Kt92t1HSYzjgX3T568G/gtocv8U7KsT2TlGUNP9U/EAduMf73SSK0V+C3ym0vFk8PluAH4H7CD5BvG3JOcXVwO/SX/un7f9Z9K+2ER61ULa3gisT9d9nZ13kY8FvkcyB8Ra4NC813wwbd8MfKDSfdFP/7yd5HD6PqAjfZzuPnolvjcB7Wn/rAcuS9vdP7v21YnsTAQ13T8uMWFmVuOqbYzAzMzKzInAzKzGORGYmdU4JwIzsxrnRGBmVuOcCKzsJJ0pKSQdleF7/HGIr78rrSbZkT7OLldslSbp1ZJ+Kik3hH28X9LX85YPlHRHke3vzK/YadXFicCycC7wc5KbaSpOiUK/6+dHREP6+H6f1+z2H9ER4IPAzRHRnd84xM80i6SYWn/+k6TqplUhJwIrq7QG0NtIboQ7J6/9xPRb+PclbZR0XV799tPTtp+ndd17asS3SPpk3j7W59d273k/Sasl3ZPWhm9O2ycpmbPgGyQVOPNv7e8v9kckXSbp58C7JZ0q6Vfpvr+XfraeOTEGFa+kuUrmCeiQ9M2eP8qS/ijpc0rmD1gj6fVp++sl3ZK23yvprZIWK51/Id3mc5I+WuCjnE9SlbWn338i6Xrg12nbrZLuVjJfwfy8/X1A0oOSfpr+G+abBaxIjwx+ln6O9ZL+Kl2/nOQLgFUhJwIrt3cBKyPiQeBpSdPy1k0FPk5S4/1Q4G2SxpKU450dEW8nKZEwGH8GzoyIacBJwJd7EgxwJPDtiJgaEY8WeO11eaeGXtuzvzSOO4FLgZPTfbcBn0jjvRqYA/wV8IaBApT0RuA9wNsiKQbXTfLHGmBvYE1EHAv8DPi7tP1rwE/T9mnA/cC/AfPSfe5Bkmiv6/Nee5LcyfpIXvN0krvwp6TLH4yI40jujP2opNem9XUuJ0kAp5D8G/XsMwccGREbgPNISsE3AMeS3NlNRDwDvCqvH62K1FU6ABt1ziUp8wtJUa9zSb6RA6yNiE4AJWWSJwF/BB6KiIfTbW4AXvmWWgIB/1fSCSRlhScAr0/XPRoRa4q89vyIaHtlR0n++E66+BaSP4a/SNv3BH4FHAU8HBG/SV+zrIR43wEcB6xL9/VqdhY1ewloTZ/fTfJHGJJa+e+DpJooSSnjZyU9JWlq+hnbI+KpPu91APCHPm1r8/oXkj/+Z6bP64HDSRLaXRGxPf1c3wGOSLeZQVLqG5L6OtcqKfx3a0R05O13G3AQ0DcmG+GcCKxs0m+DM4FjJAXJjHIh6VPpJv+Tt3k3ye9fodK8PbrofdQ6tsA255McRRwXETuUVJXs2e6FQX+Ina8RSS3/Xqc7JDXQf+ng/uIV8K2IuKTAa3bEzjovPX1SzDXA+0n+cF9bYP2f2LWfXukHSScCJ5NMnPKipLvytu/vc80GVgJExM/SpHsG8J+SvhQR3063G5u+v1UZnxqycjqb5FTMIRExKSLqSSosvr3IazYCh+ad+39P3rpHSE6LkJ5imlzg9fuQ1JffIekk4JChfYRXrCE5dfWX6fvvJemINN7Jkg5Lt8tPFP3Fuxo4W9Lr0nX7SxooztXAh9Ptc5LGpe23kJyvfzMFBm/TUzS59BRWIfsAz6RJ4CiSIx9IvvGfmJ4mGgO8O+8170jjIY17W0RcTXKqqufziiQ5PTLA57IRyInAyulckj9U+W4iOa9cUET8ieRqk5XpIO0TJKdBel67f3oa6cMkVWf7ug5olNRGcnSwcSgfIC+u7STfvG+QdB9JYjgqIv5Mciroh2m8+WMPBeNNz61fCtyR7utHJFMdFvMx4CRJvyY5ZXR0uq+XgJ+QlDbu7ue1d9B/8l0J1KVxLE4/F5HMytVCcvrrTtLTeZLGk4yb9Ez+cyLQIakdOAu4Mm0/jmSso2uAz2UjkKuPWsVJek1E/DH9VnkV8JuI+H+VjqsU6amWT0ZE0zC93x4kf6Tf3TNOUWCbqcAnIuK9ZXi/ucDEiPj8ANtdCSyPiNVDfU8bfh4jsJHg7yTNIxmQbSe5isj6kDSFZGD5lv6SAEBEtKeXjOaKHDWUJCKWlbjpeieB6uUjAjOzGucxAjOzGudEYGZW45wIzMxqnBOBmVmNcyIwM6tx/x8R2D8jGpZmcQAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "plt.errorbar(omega, Vratio, errVratio, fmt = 'ko', markersize = 5,\\\n", " linewidth = 1.8,\\\n", " markeredgecolor = 'b',\\\n", " markerfacecolor = (.49, 1, .63),\\\n", " capsize = 5)\n", "plt.xlabel('Angular Frequency (rad/s)')\n", "plt.ylabel('Voltage Ratio');" ] }, { "cell_type": "markdown", "id": "written-twenty", "metadata": {}, "source": [ "To do the actual fit, we will use the *curve_fit()* function from the *SciPy* module. This way of fitting is very nice because we use it for all types of fit models (linear, polynomial, linear-in-parameter fits, and nonlinear fits). It is capable of doing both unweighted and weighted fits and it will return uncertainties in the fit parameters via the covariance matrix.\n", "\n", "The first step is to define a function for the model that we will fit our data to. In this case, the model is a Lorentzian:
\n", "\n", "$V_\\mathrm{ratio}=\\frac{A}{\\sqrt{1+\\left(\\dfrac{\\omega}{\\gamma}\\right)^2\\left[1-\\left(\\dfrac{\\omega_0}{\\omega}\\right)^2\\right]^2}}$
\n", "\n", " > $\\omega_0$ is the resonance frequency
\n", " $A$ is the amplitude (height of the resonance peak)
\n", " $\\gamma$ is the width of the resonance
\n", " $\\omega$ is the independent variable along the horizontal axis (angular frequency in this case)\n", " \n", "- Here's how the function is defined. For $\\gamma$ we used \"width\" and for the independent variable, we used \"x\"." ] }, { "cell_type": "code", "execution_count": 20, "id": "swedish-latter", "metadata": {}, "outputs": [], "source": [ "def LRCFunc(x, A, width, w0):\n", " y = A/np.sqrt(1+(x/width)**2*(1-(w0/x)**2)**2)\n", " return y" ] }, { "cell_type": "markdown", "id": "immune-theta", "metadata": {}, "source": [ "Here is the actual command to execute the fit. At a minimum, *curve_fit()* requires as inputs the function that defines the model, the *x*-data, and the *y*-data. The statement below tells *curve_fit()* to return a list of the best-fit parameters and the covariance matrix which will be used to determine the error in the fit parameters.\n", "\n", "To give the fit a chance at being successful, we have to provide reasonable initial guesses for the fit parameters. We use the option \"p0\" for the list of starting parameters. The order must be the same as the order of the parameters defined in the function \"LRCFunc()\"." ] }, { "cell_type": "code", "execution_count": 21, "id": "defined-debate", "metadata": {}, "outputs": [], "source": [ "start = (0.7, 0.5e5, 2e5)\n", "a_fit, cov = curve_fit(LRCFunc, omega, Vratio, p0 = start)" ] }, { "cell_type": "markdown", "id": "detected-seating", "metadata": {}, "source": [ "- Here is the output of the *curve_fit()* function." ] }, { "cell_type": "code", "execution_count": 22, "id": "straight-cowboy", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "The best-fit parameters are:\n", " A = 0.7217437152596441 \n", " width = 68024.56300087317 1/s\n", " w0 = 213256.0968196423 1/s\n" ] } ], "source": [ "print('The best-fit parameters are:\\n A =', a_fit[0], '\\n',\\\n", " 'width =', a_fit[1], '1/s\\n',\\\n", " 'w0 =', a_fit[2], '1/s')" ] }, { "cell_type": "markdown", "id": "internal-construction", "metadata": {}, "source": [ "The uncertainties of the best-fit parameters are determined from the square roots of the diagonal elements of the covariance matrix. \n", "\n", "We can select the diagonal elements using:\n", "> np.diag()[0] for the (0,0) element
\n", "np.diag()[1] for the (1,1) element
\n", "np.diag()[2] for the (2,2) element\n", "\n", "Note the use of the unicode character \\u0394 to print $\\Delta$." ] }, { "cell_type": "code", "execution_count": 23, "id": "another-offset", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "The errors in the parameters are:\n", " ΔA = 0.0037048186479549886 \n", " Δwidth = 952.8034010325724 1/s\n", " Δw0 = 377.8379936633329 1/s\n" ] } ], "source": [ "print('The errors in the parameters are:\\n \\u0394A =', np.sqrt(np.diag(cov))[0],\\\n", " '\\n', '\\u0394width =', np.sqrt(np.diag(cov))[1], '1/s'\\\n", " '\\n', '\\u0394w0 =', np.sqrt(np.diag(cov))[2], '1/s')" ] }, { "cell_type": "markdown", "id": "original-delicious", "metadata": {}, "source": [ "We can now report:\n", "\n", ">$A = 0.722\\pm 0.004$
\n", "$\\gamma=68000\\pm 1000~\\mathrm{s}^{-1}$
\n", "$\\omega_0=213300\\pm 400~\\mathrm{s}^{-1}$\n", "\n", "- Next, use the function that we defined to generate the best-fit curve and plot it on top of the data." ] }, { "cell_type": "code", "execution_count": 24, "id": "adjacent-private", "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEGCAYAAABo25JHAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAA6dElEQVR4nO3deZyNdf/48dfbzNiXbG227FvJMkiLJDIYRMttGt1Uv9ykVcuNVCN1uymJEqmoOzUoS0zMaFDyzTYY+xrJliXKEsaM9++Pc9ExxsxgjmvOOe/n4zGPOee6Puc67/PhMe9zfT7X9f6IqmKMMSZ45XE7AGOMMe6yRGCMMUHOEoExxgQ5SwTGGBPkLBEYY0yQC3U7gItVqlQpveGGG9wOwxhj/MqyZcsOqGrpjPb5XSK44YYbSEpKcjsMY4zxKyKy/UL7bGjIGGOCnCUCY4wJcpYIjDEmyFkiMMaYIGeJwBhjgpwlAhOUYmJiEBHnJwSRSERecX6HnN3Xr18/t0M1xucsEZigFBMTg6qSmqpEtDtJyepjaPLiA5So9iGhBb/Fc2V1WwYNCjkvOcTExLgcvTE5y+/uIzAmJ82aBet2HeNfq+cQEqbc9aYwqtbdFDj1MyfzhVKtwwE2ftMD/eMX1iZX5brrrnY7ZGNynJ0RmKC2YgWUu2cHIWGedTlCwpRrb05BixTjX2vm0GJIMj3WzEWKVaRZs8EcP37c5YiNyXmWCExQq1cPdswuR9opASDtlLD7p2up1Pbc5FCt/X42bapLwYL/saEiE3B8OjQkIhHAcCAE+FhV/5tu/4tAtFcsNYHSqnrQl3EZExMTw4ABA4A8hBb8ltG1m1P93gNsnFaKIwfXsn7qtdw5UAgJU06dyMOa8dUpWe1aqnWow47Zz3Fj+aLETQ0lJMTtT2LM5RNfLVUpIiHAJqAlsBNYCkSp6roLtG8HPKeqzTM7bnh4uFqtIZNTVJXw8Mbs2HEjNWtGMX/+cCCB0IIzKFrmJqrfe4DV48uSt4jQY823hIQpaaeE8Y3bM/L1YkRGuv0JjMkeEVmmquEZ7fPlGUEjYIuqbnWCmAB0ADJMBEAUEOvDeIw5z9SpU1m+fCnjxj1Bt24t8XxvgbQ0z0RycnIZKjWHX8quOWeoqOw9v5KcfJMlAhMQfDlHUAbY4fV8p7PtPCJSEIgAJl9gf3cRSRKRpP379+d4oCa4pKVBXBy8/rrywgvzqFq1Bl26dDmnTUgIREZC//7QuTPs+O7ceYSNU0tRu/YpN8I3Jsf5MhFIBtsuNA7VDvi/C80NqOoYVQ1X1fDSpTMsp21Mpv6+gSyE/EXj6fr8LuKPruLPsD5s2zWMsLB85038nnlNu3Yh7N6wkNG1m5P40s2Mrt2cP3eupFOn/DZZbAKCL+cImgAxqtrKed4XQFUHZdB2KvCVqn6Z1XFtjsBcjrg46PXan3RZNP3seP/njdvzQRbj/X8PFcHNNysDB97KgQN72bRpE6GhdjuOyf0ymyPw5RnBUqCqiFQUkbxAZ2B6BsEVA+4EvvFhLMYAGd83UL7VDpKTM3+d91BRu3ZC//592bZtG1999ZXvgzbGx3yWCFQ1FXgSSADWA5NUda2I9BCRHl5NOwKzVfWYr2Ix5oyM7hvYkVCeunUv7jiRkZFUqVKF999/P+eDNOYK8+kNZao6U1WrqWplVX3T2TZaVUd7tflUVTv7Mg5jMhvv373xJ9q1C7mo8f48efLQq1cvfvrpJ1asWOHb4I3xMbuz2ASFM0XmVNOYl1CMg5v/xbVb0vjsnTKcOByBahqqmmUi8K5a+txzzwNtqV9/it1tbPyazyaLfcUmi83l6tmzJ5999hl79+6lSJEil3SMtDSI7JjK0o17qdbhALsSK1CrbGG729jkWm7dUGZMrnPy5EkmTpxIx44dLzkJgFfV0jVznauPVjG+cQdmzSpqN5kZv2NDQyaoJCYmcujQIaKjo7NunImMrj4q1+rXLK8+MiY3skRggsqUKVMoWrQoLVq0uKzjZHT10fZZZS/66iNjcgNLBCZopKamMn36dCIjI8mbN+8lHSOzq49+27zooq8+MiY3sDkCEzQWLFjAgQMH6NSp0yUfIyYm5uwfee/CdFryTX7+4312/XGCsLCwHIrYmCvDzghM0JgyZQr58+cnIiIiR47nfbdxv3512L//N2bNmpUjxzbmSrJEYILC6dOnmTJlChERERQqVCjHjx8REUHp0qUZP358jh/bGF+zRGCCQlJSErt27aJjx44+OX5YWBgPPvggM2bM4PDhwz55D2N8xRKBCQozZ85ERGjbtq3P3iMqKooTJ07wzTdWP9H4F0sEJigkJCTQqFEjSpYs6bP3aNSoCaVLP8LrrytxcZ7JZGP8gSUCE/AOHTrEkiVLaNWqVY4f23vBm4JXzSat+EBKd7yZrs/vIn/R+LP1h+xyUpObWSIwAS8xMZHTp0/7LBGoKjNmpHF9jSb0WDOXFkNW0mPNXK6vfiszZmSvmJ0xbrJEYAJefHw8xYoVo1GjRj57Dys5YfyZJQIT0FSVhIQEWrRo4dMlJTMqOfFrvJWcMP7BEoEJaOvWrWPXrl05dhNZepmVnNizyUpOGP9gJSZMQEtISADwyfwAXLjkxAF6U7bK96xcaZcOmdzPEoEJaAkJCdSsWZNy5cr5/L3OlJyIjIQiRSrw7LMrWL9+PTVr1vT5extzOWxoyASsv/76ix9++MFnZwOZeeCBBxARJkyYcMXf25iL5dNEICIRIrJRRLaISJ8LtGkmIskislZEfvBlPCa4zJ8/n5MnT7qSCK6//nqaNWtGbGws/rYcrAk+PksEIhICjARaA7WAKBGpla7NVcAHQHtVrQ084Kt4TPBJSEggf/783Hnnna68f1RUFJs3b2b58uWuvL8x2eXLM4JGwBZV3aqqKcAEoEO6Ng8BU1T1VwBV3efDeEyQSUhIoGnTphQoUMCV97/vvvsICwsjNjbWlfc3Jrt8mQjKADu8nu90tnmrBhQXke9FZJmI/DOjA4lIdxFJEpGk/fv3+yhcE0h+/fVX1q9f78qw0BklSpSgVatWTJw4kdOnT7sWhzFZ8WUikAy2pR8sDQUaAG2BVsArIlLtvBepjlHVcFUNL126dM5HagKOry8bza6oqCh27tzJggULXI3DmMz4MhHsBLyv2SsL7M6gTbyqHlPVA8B84GYfxmSCREJCAmXKlKFWrVpZN/ah9u3bU6BAARseMrmaLxPBUqCqiFQUkbxAZ2B6ujbfAHeISKiIFAQaA+t9GJMJAqmpqSQmJhIREYFIRiemV07hwoVp3749X3/9NadOnXI1FmMuxGeJQFVTgSeBBDx/3Cep6loR6SEiPZw264F4YBWwBPhYVdf4KiYTHJYsWcKff/7p+rDQGVFRURw4cIA5c+a4HYoxGfLpncWqOhOYmW7b6HTP3wLe8mUcJjicKfEwdOgJRNpx110t3A4J8KxnXKxYMWJjY31W88iYy2F3Fhu/5r0wTP6i8XR9fhcnG5akeNVRXFdhca5YGCZfvnx06tSJqVOncvz4cdfiMOZCLBEYv+YvC8NERUVx5MgRZs6cmXVjY64wSwQmIOT2hWHuuusurr76art6yORKlghMQMhoYZgdCeVzzcIwoaGhPPjgg8TFxXH48GG3wzHmHJYIjF/LbGGY3Rt/ylULw0RFRXHy5EmmTZvmdijGnEP8rTJieHi4JiUluR2GyYWSk1dTr15f7r13AI891oDWrT1rBOQWqkrFihWpVauWzRWYK05ElqlqeEb77IzABIzExATgW9577xoiI3NXEgAQETp37szs2XP48ss/GTgQ4uI8l70a4yZLBCZgJCQkULt2bcqWLet2KOc5M4Q1ePBbSL5veGrAUWYdWUnX53eRv2h8rrjM1QQvSwQmIBw7doz58+fnmruJ0ztzmev06alcVf7mXHuZqwlOlghMQPjhhx9ISUnJtYngjORkoVr7A7n2MlcTnCwRmIBwZjWyO+64w+1QMuW5zLVsrr3M1QQnSwQmICQkJNCsWTPXViPLivdlrns2Lc7Vl7ma4GOJwPi97du3s3Hjxlw9LHRmjkA1jROHI3ig+QoWvvUVMU+d5MThCFRtjsC4xxKB8Xvx8fEAflPZMyQE3nyzCWFhQ9i69b1cd5mrCT6WCIzfi4+Pp0KFClSvXt3tULKtZMmStG/fnvHjx5OSkuJ2OCbIWSIwfi0lJYU5c+bkitXILtYjjzzCgQMHmDVrltuhmCBnicD4tYULF3LkyBG/GRby1qpVK6655hrGjRvndigmyFkiMH4tPj6e0NBQ7r77brdDuWihoaE8/PDDfPvtt+zbt8/tcEwQs0Rg/Fp8fDy33347RYoUcTuUS9KtWzdSU1P58ssv3Q7FBDFLBMZv7dmzh+TkZL8cFjqjdu3aNGzYkE8//dTtUEwQ82kiEJEIEdkoIltEpE8G+5uJyJ8ikuz8vOrLeExgmT17NuA/l41eSLdu3Vi5ciXLly93OxQTpHyWCEQkBBgJtAZqAVEiUiuDpj+qal3n53VfxWMCT3x8PNdeey116tRxO5TL8tBDD1GgQAE+/PBDt0MxQcqXZwSNgC2qulVVU4AJQAcfvp8JImlpacyePdsvLxtN76qrrqJz58588cUXtoylcYUvE0EZYIfX853OtvSaiMhKEZklIrUzOpCIdBeRJBFJ2r9/vy9iNX5m6dKlHDx40O+Hhc7o0aMHx44ds0lj4wpfJoKMvqalXxdzOVBBVW8G3gOmZXQgVR2jquGqGl66dOmcjdL4pfj4ePLkyUOLFi3cDiVHNGzYkLp16zJ69Gj8bflY4/98mQh2AuW8npcFdns3UNXDqnrUeTwTCBORUj6MyQSI+Ph4GjVqRMmSJd0OJUeICD169GDlypUsXrzY7XBMkPFlIlgKVBWRiiKSF+gMTPduICLXijPAKyKNnHh+92FMJgDs3buXJUuW0KZNG7dDyVEPPfQQhQsXtkljc8X5LBGoairwJJAArAcmqepaEekhIj2cZvcDa0RkJTAC6Kx2XmyyMHPmTFSVdu3auR1KjipSpAjR0dFMmDCBQ4cOuR2OCSKSnb+7ItIeaOo8/UFVZ/g0qkyEh4drUlKSW29vcoFOnTqRlJTE9u3b/f6KofSSk5OpV68e77zzDs8995zb4ZgAIiLLVDU8o31ZnhGIyCDgGWCd8/O0s82YK+7EiRPMnj2byMjIgEsCAHXr1uX222/nvffeIy0tze1wTJDIztBQW6Clqo5V1bFAhLPNmCvu+++/59ixY0RGRrodis88++yzbNu2jRkzXDvxNkEmu3MEV3k9LuaDOIzJlri4OAoWLEjz5s3dDsVnOnToQIUKFXj33XfdDsUEiewkgkHAChH5VEQ+A5YB//FtWMacT1WZMWMGLVu2JH/+/G6H4zOhoaE8/fTT/PDDD6xYscLtcEwQyDIRqGoscAswxflpoqoTfB2YMemtXr2aX3/9NaCHhc547LHHKFy4MMOHD3c7FBMELpgIRKSG87s+cB2eG8R2ANc724y5YtLS4K231gP9yZu3E4E+j1qsWDEeeeQRYmNj+e2339wOxwS4zM4Ieju/h2bw87aP4zLmrLQ0iOyYSkJSU2596QFeGR5CZMfUgE8GTz31FKdOnWLkyJFuh2IC3AUTgap2dx62VtW7vH+AwLql0+RKMTExiAihoZEs2biXx1clcvfgZLosms6SDXsJDfVcQhoTE+N2qD5RtWpV7r33Xt5//32OHDnidjgmgGVnsvinbG4zJkfFxMSgqrz+ehw17j1ISJjn5seQMKXmfYcYODAOVQ3YRADQt29f/vjjD0aPHu12KCaAZTZHcK2INAAKiEg9Eanv/DQDCl6pAI2pVw+2zChN2inPDWRpp4QdCeWpW9fduK6Ehg0b0rJlS4YOHcqJEyfcDscEqMzOCFrhmQsoC7zD3/MDvYF+vg/NBLszQ0Pt2oVwcHsyo2s3J/Glmxlduzm7N/5Eu3YhAT00dEbfvn3Zu3cv48aNczsUE6CyrDUkIvep6uQrFE+WrNZQ8Pn000955JHHeOedDRw7VpW6daF1awgJcTuyK0NVufXWW/ntt9/YvHkzoaGhbodk/FBmtYay/B+lqpNFpC1QG8jvtd3WFzZXxOTJkylfvizPPluFACwvlCURoV+/frRv357Y2Fgefvhht0MyASY7RedGA/8AnsKz6tgDQAUfx2UMAIcPH2b27Nl06tQpIIvMZVfbtm2pU6cOAwcOJDU11e1wTIDJzlVDt6rqP4FDqjoAaMK5K48Z4zMzZ84kJSWF++67z+1QXJUnTx5ef/11Nm/ezOeff+52OCbAZCcRHHd+/yUi1wOngIq+C8mYv02ePJlrrrmGJk2auB2K69q3b094eDgDBgzg5MmTbodjAkh2EkGciFwFvIVnsflfAKs1ZHzur7/+YubMmXTs2JGQYJkZzoSI8MYbb7B9+3Y++eQT0tIgLg4GDvT8DvQ7rY3vZGuFsrONRfLhmTBOVdVjPosqE3bVUPCYOHEinTt3Zu7cudx1111uh5MrqCp33nknmzdv5aYGP7NxzwnK3bODHbPLUatMIeKmhgbN1VTm4lzyCmUiUkZEwp3F58GzFsG/gc05HKMx54mNjeW6666jadOmWTcOEgMGDODHH3/kt9/qsmzzAbosmk7zQcuDpuyG8Y3M7ix+FkgG3gMWiUhXPIvQFwAaXIngTPD6448/mDVrFv/4xz9sWMjLmbIbVao8SPUOvwdl2Q2T8zI7I+gOVFfVJsC9wEdAW1V9TlX3ZOfgIhIhIhtFZIuI9MmkXUMRSROR+y8meBO4pkyZQkpKCg899JDboeRKzzzTlI3flAzKshsm52WWCE6o6kEAVf0V2KSqi7J7YBEJAUYCrYFaQJSI1LpAu8FAwsUEbgJbbGwslStXJjw8wyHNoHWm7MZTT1Xm8M7VQVt2w+SszO4sLisiI7yeX+39XFWfzuLYjYAtqroVQEQmAB2AdenaPQVMBhpmO2oT0H777Tfmzp1Lv379gvomsozExMSc/SO/e/deKld+kj/i7uWzd6Jp3boMISF26ZC5eJklghfTPV92kccug2dFszN2Ao29G4hIGaAj0JxMEoGIdMczVEX58uUvMgzjbyZNmsTp06eJiopyO5Rc7frrr+HVV+vTr18XChcuQ0hIM7dDMn7qoi4fvagDizwAtFLV/+c8fxhopKpPebX5ChiqqotE5FMgTlW/zuy4dvlo4GvSpAl//fUXK1eudDuUXO/48ePUqFGDq666imXLlllBOnNBl3z56GXaybmlKMoCu9O1CQcmiMgvwP3AByJyrw9jMrncpk2bWLRoEdHR0W6H4hcKFCjAsGHDWLVqFe+//77b4Rg/5ctEsBSoKiIVnfsQOgPTvRuoakVVvUFVbwC+Bp5Q1Wk+jMnkcp9++il58uSxCpsXoWPHjrRu3ZpXXnmFXbt2uR2O8UM+SwSqmgo8iedqoPXAJFVdKyI9RKSHr97X+K+0tDT+97//ERERwXXXXed2OH5DRHj//fdJTU2ld+/ebodj/FB2ylBXE5E5IrLGeV5HRPpn5+CqOlNVq6lqZVV909k2WlXPW4BVVbtlNT9gAltiYiK7du3ikUcecTsUv1OpUiX69evHpEmTmD17ttvhGD+TnTOCj4C+eKqOoqqr8AzzGJOjxo0bR4kSJWjXrp3bofill156iapVq9KzZ0+OHXOlFJjxU9lJBAVVdUm6bbYyhslRhw4dYtq0aURHR5MvXz63w/FL+fLl46OPPmLr1q28/PLLbodj/Eh2EsEBEakMKIBTBiJbJSaMya7Y2FhOnjxpw0KX6c477+TJJ59kxIgR/Pjjj26HY/xEdhavrwSMAW4FDgHbgC6q+ovPo8uA3UcQmMLDwzl16hTJycl2N/FlOnbsGHXq1EFEWLVqFQULFnQ7JJMLXNZ9BKq6VVVbAKWBGqp6u1tJwASmpUuXsmzZMrp3725JIAcUKlSITz75hJ9//tmGiEy2ZHkbooj0Tvcc4E9gmaom+yYsE0xGjRpFoUKF7N6BHNSsWTN69erFu+++S9u2bWnRooXbIZlcLDtzBOFADzy1g8rgqfnTDPhIRF7yXWgmGBw6dIgJEyYQHR1N0aJF3Q4noAwZMoSaNWvy8MMPs3//frfDMblYdhJBSaC+qj6vqs/jSQylgaZANx/GZoLAZ599xvHjx+nZs6fboQScggULMmHCBA4dOkS3bt3wVV0x4/+ykwjKAylez08BFVT1OHDSJ1GZoKCqjB49mltuuYW6tqKKT9SpU4e3336bmTNnMmLEiKxfYIJSdkoVfolnqcpvnOftgFgRKcT5awsYk23z5s1j48aNfPbZZ26HEtB69erF7Nmzeemll7j99ttp0MBWmjXnylYZahEJB24DBFigqq5dv2mXjwaOTp068cMPP7Bz504KFCjgdjgB7cCBA9SvX588efKQlJREqVKl3A7JXGGXXYba+cMfC0wB9omIrQ5jLllaGnz00R6mTr2Ru+8eRt68lgR8rVSpUkyePJnffvuNqKgoUlOtOID5W3aKzrUXkc14biT7wfk9y9eBmcByZq1dkRDyF42nz9unafLifcxZeTf5i8YjYmvt+lrDhg354IMPSExMpH//bNWNNEEiO3cWr8SzlGSiqtYTkbuAKFXtfiUCTM+GhvxbXBz0evVPuiyeTkiYknZKGN+4AyNfL0pkpNvRBYcePXrw4YcfMmnSJB544AG3wzFXyOUODZ1S1d+BPCKSR1XnAXVzMkATPFasgLItfyUkzPMFJCRMKdfqV5KT3Y0rmAwfPpwmTZrQtWtXlixJX0/SBKPsJII/RKQwMB/4QkSGY9VHzUU6MzT06quRbJhWgrRTnlISaaeE9ZOL88orkTY0dIXky5ePadOmce2119KuXTu2bNlGXBwMHOg5Y0tLcztCc6VlZ2ioEHAcT9KIBooB41X1oO/DO58NDfm3zz//kkd7lODqyo2o2HY3OxLKU6tsQeKmhhIS4nZ0wSMmJoYBAwYAeQgt+C1Fy95E9Q4H2PhNKQ7vXE3qX22B07z22muWnANEZkND2UkEg1X131ltu1IsEfiv06dPU6dOHU6fFv7735WsWpWHunWhdWssCbhk0KDVvP1pCXqsmWtzNgHucucIWmawrfXlhWSC0TfffMPatWvp378v7dvnoX9/iIy0JOCm1NSbqNHxoM3ZBLkLJgIR6Skiq4HqIrLK62cbsOrKhWgCgaryxhtvUKVKFR588EG3wwl658zZTLU5m2B3waEhESkGFAcGAX28dh3J7vyAiEQAw4EQ4GNV/W+6/R2AgcBpPBPQz6rqgsyOaUND/mnWrFm0adOGTz75hEcffdTtcIwjLQ0iO6aybudflG25nQ1TS1Aszx9sWlOL0FBbGyKQXNIcgYiUyOygWSUDEQkBNuEZWtoJLMVz/8E6rzaFgWOqqiJSB5ikqjUyO64lAv+jqtx2223s2rWLzZs3kzdvXrdDMl7S0mDWLFixQklK+pjp03vQr18f3njjDVsoKIBklggyKzq3DGedYjw1hrwpUCmL920EbFHVrU4QE4AOeBWqU9WjXu0Leb2fCSDz5s1j4cKFjBw50pJALhQS4pmriYwUVP8f//rXUv7zn/+QlpbGoEGDLBkEgQsmAlWteJnHLgPs8Hq+E2icvpGIdMQz/HQ10DajA4lIdzwL4lC+vJU58ieqSr9+/ShbtqwNCfkBEWH06NGEhoYyePBgUlJSGDp0qCWDAJedMtSISHs8C9EAfK+qcdl5WQbbzvvGr6pTgaki0hTPfMF5a+qp6hhgDHiGhrITs8kdpk+fzuLFi/noo4/Inz+/2+GYbMiTJ8/Zs7dhw4aRkpLCiBEjyJMnWzUqjR/KzprF/wUaAl84m54RkdtUtW8WL90JlPN6XhbYfaHGqjpfRCqLSClVPZBVXCb3S0tL4+WXX6ZatWp069bN7XDMRRARhg0bRlhYGG+//TaHDx/mk08+ISwszO3QjA9k54ygDVBXVU8DiMhnwAogq0SwFKgqIhWBXUBn4CHvBiJSBfjZmSyuD+QFfr+4j2Byq9jYWNauXcvEiRMJDc3WyafJRUSEIUOGUKxYMV555RX27dvH119/TeHChd0OzeSw7J7rXeX1uFh2XqCqqcCTQAKwHs8VQWtFpIeI9HCa3QesEZFkYCTwD7WFVQNCSkoKr776KvXq1eP+++93OxxziUSE/v378/HHH5OYmEizZs3Yu3ev22GZHJadr2mDgBUiMg/PuH9Tsj4bAEBVZwIz020b7fV4MDA429EavzF69Gi2bdvGrFmzbGw5ADz22GNce+21PPjgg9x2223ExcVRo0amV3obP5LZfQTvA1+q6k8ich2eeQIBFqvqb1cwxnPYfQS53++//07VqlUJDw8nISHBrjgJIIsXL6Z9+/acOHGCCRMm0Lq1VZvxF5daa2gzMFREfgGeBX5V1W/cTALGP8TExPDnn3/yzjvvWBIIMI0bN2bp0qVUqlSJtm3b8vbbb2Ojuf7vgolAVYerahPgTuAgME5E1ovIqyJS7YpFaPzK2rVrGTVqFD169ODGG290OxzjA+XLl2fBggXcf//9vPjii3Tt2pXjx4+7HZa5DFkO3qrqdlUdrKr18Fz10xHP5K8x51BVevfuTZEiRZxa9yZQFSpUiIkTJzJw4EA+//xzGjduzMaNG90Oy1yi7CxeHyYi7UTkCzyL1m/Cc7WPMeeYMWMGs2fP5rXXXqNUqVJuh2N87MwVRbNmzWLPnj00aNCAL7/80u2wzCXIbLK4JRCFp+zDEmACME1Vj1258M5nk8W509GjR6lVqxbFihVj+fLlduNRkNm5cydRUVEsWLCAxx9/nOHDh1OgQAG3wzJeLnWyuB+wEKipqu1U9Qu3k4DJvWJiYtixYwcffvihJYEgVLZsWebNm0efPn346KOPqF+/PklJSaSlYesh+4HMJovvUtWP3Fqb2PiP5ORk3n33Xbp3786tt97qdjjGJaGhoQwaNIguXbqwYcMGGjZsTP6i8XR9fhezjqyk6/O7yF80HpEQW/Qml8lyzeLcxoaGcpe0tDRuvfVWfvnlFzZs2EDx4sXdDsnkAocOHaJTp7Gs2t3Z1kPOJS53zWJjLmjkyJEsWbKEYcOGWRIwZxUvXpzmzZ8/fz3ke2w95NzIEoG5ZJs3b6ZPnz60bt2aqKgot8MxuUSm6yFPsfWQcyMbGjKXJC0tjTvuuIMNGzawZs0arr/+erdDMrmM93rI5Vr9yi8zr+fwztUcOdicxx57hMGDB1OyZEm3wwwaNjRkctzQoUNZuHAh7733niUBk6GQEIibGsrI14sSUehGRr9Zgp3bwnnxxef59NNPqVq1Ku+99x6nTp1yO9SgZ2cE5qKtXbuW+vXrExkZyddff231hMxFW716Nc899xxz5syhZs2aDBs2jFatWrkdVkCzMwKTY06cOEF0dDTFihVj1KhRlgTMJbnpppv47rvvmDZtGikpKURERBAZGcnatWvdDi0oWSIwF+WFF15g5cqVjBs3jquvvtrtcIwfExE6dOjA2rVrGTJkCPPnz+emm26ia9eubNu2ze3wgoolApNtU6ZMYeTIkfTu3Zu2bdu6HY4JEPny5ePFF19k69atPP/880yaNInq1avTq1cv9uzZ43Z4QcHmCEy2bN++nbp161K1alUWLFhA3rx53Q7JBKhdu3YxcOBAPvnkE8LCwujZsyfPP/+8XZRwmWyOwFyWlJQUoqKiOH36NBMmTLAkYHyqTJkyjB49mvXr13Pffffx7rvvUrFiRXr27GlDRj5iicBk6ZlnnmHhwoV8/PHHVKpUye1wTJCoUqUKn3/+OZs2baJbt26MHTuWqlWr8s9//pN169a5HV5A8WkiEJEIEdkoIltEpE8G+6NFZJXz85OI3OzLeMzF+/jjjxk9ejT//ve/eeCBB9wOxwShypUr8+GHH7J161aefvppJk+eTO3atYmIiCA+Pp7Tp0+7HaL/U1Wf/AAhwM9AJSAvsBKola7NrUBx53FrYHFWx23QoIGaK2PhwoWaN29eveeeezQ1NdXtcIxRVdX9+/frwIED9dprr1VAa9SooaNGjdKjR4+6HVquBiTpBf6u+vKMoBGwRVW3qmoKnoVtOqRLQj+p6iHn6SKgrA/jMRdhz549dOrUibJlyxIbG0tISIjbIRkDQKlSpejfvz/bt2/n888/p3DhwvTs2ZNy5crxwgsvsGnTJrdD9Du+TARlgB1ez3c62y7kMTxLYZ5HRLqLSJKIJO3fvz8HQzQZOXr0KJGRkRw+fJipU6dSokQJt0My5jx58+alS5cuLFmyhAULFtCiRQuGDx9O9erVadasGV9++SUnTpxwO0y/4MtEkNEtpxleqyoid+FJBP/OaL+qjlHVcFUNL126dA6GaNJLTU2lc+fOrFixihde+J5vvqljK0uZXE1EuO2225g0aRI7duxg0KBB7Nixg+joaMqUKcNzzz13zh3Ltmra+Xx2H4GINAFiVLWV87wvgKoOSteuDjAVaK2qWZ7T2X0EvqOq9OrVi1GjPuTGej9zWIpT7p4d7JhdjlplChE3NRQbITL+4PTp08ybN48xY8YwdepUTp06Rf369XnooYf59rvH+Xl/atD933brPoKlQFURqSgieYHOwPR0gZUHpgAPZycJGN+65557GDVqFNCa3cfC6LJoOs0HLafLouks2bCX0FCrI2/8Q548ebj77ruZOHEiO3fupFWrVixfvpwXXkhk5bY/7P92OqG+OrCqporIk0ACniuIxqrqWhHp4ewfDbwKlAQ+cIqXpV4oYxnfGjt2LImJiTz44IPUrh3L7OPrzllZquZ9h4goFEf//i4HasxFuvrqq4mPjwfg6af3s6zAb+f8367R8XfK//pfPv54IoUKFXIzVNf49D4CVZ2pqtVUtbKqvulsG+0kAVT1/6lqcVWt6/xYEnDBV199xeOPP07lypWZNGkSr73WnvVTip+7stRkW1nK+K8zq6a9994jbJh27qppG6aWZMKEPhQuXJhatWoRGxvL4cOHXY74yrJaQ0EuPj6e9u3b06hRIxISEihUqNB5K0vtSChPrbIFg2Ic1QS2jP5v1yxTgBef/T+mTv2aKVOmsGfPHvLmzUurVq3o2LEjbdq04ZprrnE79MuW2RyBJYIgNn/+fCIiIqhZsyZz586lWLFiZ/elpcGsWZCcDHXrQuvWWBIwASGz/9unT59m4cKFfP3110yePJkdOzxXwDds2JDIyEgiIyOpV6+eX67DYYnAnOf777+nbdu2VKhQgR9++AG7LNeYc6kqK1euJC4ujm+//ZbFixejqlx//fW0adOGyMhImjdvTpEiRdwONVssEZhzfPfdd3To0IFKlSoxZ86cgDjtNcbX9u3bx6xZs4iLi2P27NkcPnyY0NBQGjduTIsWLWjZsiWNGjUiLCzM7VAzZInAnDVz5kw6depEjRo1+O677+xMwJhLkJKSwv/93//x3XffkZiYSFJSEqpKkSJFaNasGS1atKBFixbUrFkz1wwjWSIwAEydOpV//OMf1KlTh9mzZ1vpCGNyyMGDB5k3bx6JiYl89913/PzzzwBcc8013HHHHTRt2pQ77riDm266ybW6XZYIDB9++CFPPPEEDRs2JD4+nquuusrtkIwJWNu2bSMxMZH58+czf/58fv31VwCKFSvGbbfdRtOmTWnatCkNGjS4Ygs9WSIIYqrKgAEDGDBgAG3atGHSpElBe9OMMW7Zvn07P/74I/Pnz+fHH39kw4YNABQoUIDw8HAaN27MLbfcQuPGjSlb1jdFmC0RBKm0tDSeeOIJxowZQ7du3RgzZkyuncgyJpjs27ePBQsWMH/+fBYtWsSKFStISUkBPEt1eieGBg0a5MiXN0sEQejIkSNER0czY8YM+vXrxxtvvJFrJq2MMec6efIkycnJLF68mEWLFrF48WK2bt0KQEhICLVr1+Hqq7uSL98t3HlnUf71r/IULXpxycESQZD55ZdfaN++PevWrWPEiBE88cQTbodkjLlI+/fvZ/HixSxcuIRP/teO1EJlqNZ+Pxu/KcWRnWuoXP55GjSoS7169ahfvz716tXLdO4vs0Tgs6Jzxh0LFiygY8eOpKamEh8fT4sWLdwOyRhzCUaOHMmAAQOAtpSo9i96rJ5DSJhy15vC6NrN2bDhBjZs+IIvvvji7GsqVqxI3bp1qVOnDjfddBN16tShUqVKWV6pZGcEAWTcuHH861//omLFisyYMYNq1aq5HZIx5jINHAgJf62h+aDlZ7fN7VufiEI30r+/Z75hxYoVLF++nOXLl7Nq1Sq2bNnC6dOnAShYsCA33ngjS5YscWU9AnOFnDhxgh49evDoo4/SrFkzFi1aZEnAGD93pmLqq69GZloN+IMPPqBVq1b07duXr776io0bN3LkyBGWLl3KJ598Qvfu3SlcuHCm72VnBH5u27Zt3H///Sxfvpw+ffowcOBAQkNtxM+YQJFT1YBtjiBATZ06nYcf/pK0tE688soHvPZaY6sQakyACQmBuKmhzJpVlOTkG6n7es5XA7YzAj+UkpJCv36vMHzUXRQvfzPVOvzOju+CZ+1VY8zFc2vNYuMDvXr1Il++fAwdupaiZW/i8VWJNP+vrb1qjLl0NjTkJ1SVUaNGMW7cOEqVKkXLlm/xa4VDtq6wMeay2RmBH9i3bx/t2rWjV69elClThgMHDhAb+6KtK2yMyRE+TQQiEiEiG0Vki4j0yWB/DRFZKCInReQFX8bij1SV2NhYatWqRWJiIiNGjGDTpk2oKqmpcTSqfg3jG3dgbt/6jG/cgUY1riE1NQ5VtURgjMk2nw0NiUgIMBJoCewElorIdFVd59XsIPA0cK+v4vBXu3fvpmfPnkyfPp3GjRszduxYatWqdXb/lbiSwBgTHHw5R9AI2KKqWwFEZALQATibCFR1H7BPRNr6MA6/oqqMGzeO3r17c/LkSYYOHcozzzyT4S3iISEQGen5McaYS+XLoaEywA6v5zudbRdNRLqLSJKIJO3fvz9HgsuNNm/ezD333MNjjz1G3bp1Wb16Nb1793ZtRSNjTHDwZSLIqObxJd20oKpjVDVcVcMDcY3d48eP8+qrr56pB8LIkSOZO3cuVapUcTs0Y0wQ8OXQ0E6gnNfzssBuH76fX/r222956qmn2LZtG9HR0bz11ltcd911bodljAkivjwjWApUFZGKIpIX6AxM9+H7+ZXt27fTsWNHIiMjyZ8/P3PnzmX8+PGWBIwxV5zPzghUNVVEngQSgBBgrKquFZEezv7RInItkAQUBU6LyLNALVU97Ku43HbkyBEGDx7M0KFDyZMnD4MHD+bZZ5+9YgtYG2NMej69s1hVZwIz020b7fX4NzxDRgEvLS2NcePG0b9/f/bu3Ut0dDT/+c9/KF++vNuhGWOCnJWY8LG0NBg0KJn33lvAvn3x3HprNaZPn06jRo3cDs0YYwBLBD6VnLya1u1PklLgOmp0bUpYQhRFShalQYMwt0MzxpizrNaQD2zZsoUuXbpQr14/ThW4nh5r5nL3kGS6Js1k/c7jzJrldoTGGPM3SwQ5aNeuXfTo0YPq1as7C0rXpVqH/edUCC3TYhvt2vW3wnDGmFzDEkEOOHDgAC+88AJVqlRh7NixNGjQwNmzgo3flDqnQujGaaWAZLdCNcaY89gcwWX4/fffGTZsGCNGjODo0aM8/PDDxMTEULFiReDvtUbHN+5wdq3RRjUKErc+zorDGWNyDUsEl2Dfvn0MHTqUkSNH8tdff3HfffcRExND7dq1z2lnFUKNMf7AEsFF2LNnD2+99RajR4/mxIkTdO7cmZdffvm8BODNKoQaY3I7SwTZsHPnToYMGcKYMWNITU0lOjqafv36Ub16dbdDM8aYy2aJIBMbN25k6NChfPbZZ5w+fZquXbvSt29fKleu7HZoxhiTYywRZGDhwoUMGTKEb775hnz58vHoo4/y73//mxtuuMHt0IwxJsdZInCcPn2ab7/9liFDhrBgwQKKFy/Oyy+/zFNPPcXVV1/tdnjGGOMzQZ8IUlJS+PLLL3nrrbdYt24d5cuX59133+Wxxx6jcOHCbodnjDE+F7SJ4M8//+Sjjz5i2LBh7N69mzp16jB+/HgefPBBwsKsFpAxJngEfCJIS4NZs2DFCqhXD6pU2cwHH7zHuHHjOHr0KM2bN2fs2LHcc889iGS0uqYxxgS2gE4EZ+7sXbfrGOVa/sr7L13Noe1b0ZQxREU9yDPPPONVDsIYY4JTQCeCadNOsuqXI3RbNouQMOXOgcLYehEMfmk3//xnCbfDM8aYXCEgi841adIEEeH++9+gQsSuc6p/Vmyzm65d30FEaNasmbuBGmNMLhBQiWDp0qVER0eTlJSEiHDLLfnZMbvcOdU/dyVWZMaMN1BVvv/+e3cDNsaYXMDvh4ZSU1OZOnUq7777Lj/99BNFihThySef5KmnnqJChUrnVf+sVbYgrVu7HbUxxuQePk0EIhIBDAdCgI9V9b/p9ouzvw3wF9BNVZdn59gHDx7k448/5v3332fHjh1UrlyZ4cOH061bN4oWLXq2nVX/NMaYzImq+ubAIiHAJqAlsBNYCkSp6jqvNm2Ap/AkgsbAcFVtnNlxb7zxRm3atCmfffYZf/31F82bN+eZZ56hbdu2hNhfeGOMyZCILFPV8Iz2+fKMoBGwRVW3OkFMADoA67zadAD+p55stEhErhKR61R1z4UOunbtWrZs2UJ0dDTPPPMMderU8eFHMMaYwOfLRFAG2OH1fCeeb/1ZtSkDnJMIRKQ70N15evLkyZNrxo4dy9ixY3M24sBRCjjgdhC5nPVR1qyPMudv/VPhQjt8mQgyuk03/ThUdtqgqmOAMQAiknSh0xvjYX2UNeujrFkfZS6Q+seXl4/uBMp5PS8L7L6ENsYYY3zIl4lgKVBVRCqKSF6gMzA9XZvpwD/F4xbgz8zmB4wxxuQ8nw0NqWqqiDwJJOC5fHSsqq4VkR7O/tHATDxXDG3Bc/noI9k49BgfhRxIrI+yZn2UNeujzAVM//js8lFjjDH+IaBKTBhjjLl4lgiMMSbI+VUiEJEIEdkoIltEpI/b8eQ0ERkrIvtEZI3XthIi8p2IbHZ+F/fa19fpi40i0sprewMRWe3sG+GU8kBE8onIRGf7YhG5wes1XZ332CwiXa/QR74oIlJOROaJyHoRWSsizzjbrY8cIpJfRJaIyEqnjwY4262P0hGREBFZISJxzvPg7SNV9YsfPBPOPwOVgLzASqCW23Hl8GdsCtQH1nhtGwL0cR73AQY7j2s5fZAPqOj0TYizbwnQBM99GrOA1s72J4DRzuPOwETncQlgq/O7uPO4uNv9kUH/XAfUdx4XwVPCpJb10Tl9JEBh53EYsBi4xfoow77qDXwJxDnPg7aPXP/HuIh/tCZAgtfzvkBft+Pywee8gXMTwUbgOufxdcDGjD4/nquzmjhtNnhtjwI+9G7jPA7Fc1ekeLdx9n2Ipy6U6/2RRV99g6eWlfVRxv1TEFiO545+66Nz+6YsMAdozt+JIGj7yJ+Ghi5UjiLQXaPOvRXO76ud7RfqjzLO4/Tbz3mNqqYCfwIlMzlWruWcatfD843X+siLM+SRDOwDvlNV66PzvQu8BJz22ha0feRPiSBb5SiCyIX6I7N+upTX5DoiUhiYDDyrqocza5rBtoDvI1VNU9W6eL71NhKRGzNpHnR9JCKRwD5VXZbdl2SwLaD6yJ8SQbCWo9grItcBOL/3Odsv1B87ncfpt5/zGhEJBYoBBzM5Vq4jImF4ksAXqjrF2Wx9lAFV/QP4HojA+sjbbUB7EfkFmAA0F5HxBHMfuT02dRFjeqF4JlYq8vdkcW234/LB57yBc+cI3uLcCawhzuPanDuBtZW/J7CW4pkgPDOB1cbZ3otzJ7AmOY9LANvwTF4Vdx6XcLsvMugbAf4HvJtuu/XR331RGrjKeVwA+BGItD66YH814+85gqDtI9f/IS7yH60NnitFfgZedjseH3y+WDwluE/h+ebwGJ5xxTnAZud3Ca/2Lzt9sRHnagVneziwxtn3Pn/fQZ4f+ApPSY8lQCWv1zzqbN8CPOJ2X1ygf27Hcxq9Ckh2ftpYH53TR3WAFU4frQFedbZbH2XcX834OxEEbR9ZiQljjAly/jRHYIwxxgcsERhjTJCzRGCMMUHOEoExxgQ5SwTGGBPkLBGYHCciHUVERaSGD9/j6GW+/nunkmSy83N/TsXmNhEpICI/iEjIZRyjm4i87/X8OhGZnUn7RO9qnca/WCIwvhAFLMBzI43rxCOj/+vRqlrX+fk63Wsu+Y9oLvAoMEVV07w3XuZnisBTSO1CPsdTcdP4IUsEJkc5dYBuw3MzXGev7c2cb+Ffi8gGEfnCq3Z7G2fbAqem+5n68DEi8oLXMdZ413U/834iMkdEljt14Ts4228Qz7oFH+CpwOl9W/+FYv9FRF4VkQXAAyJyj4gsdI79lfPZzqyLcVHxikgX8awTkCwiH575oywiR0XkTfGsH7BIRK5xtl8jIlOd7StF5FYRGSjOGgxOmzdF5OkMPko0nsqsZ/p9noh8Cax2tk0TkWXiWa+gu9fxHhGRTSLyg/Nv6C0CmOWcGcx3PscaEbnD2T8dzxcA44csEZicdi8Qr6qbgIMiUt9rXz3gWTz13SsBt4lIfjyleFur6u14SiRcjBNAR1WtD9wFDD2TYIDqwP9UtZ6qbs/gtV94DQ2VPHM8J45EoD/Qwjl2EtDbifcjoB1wB3BtVgGKSE3gH8Bt6ikGl4bnjzVAIWCRqt4MzAced7aPAH5wttcH1gKfAF2dY+bBk2i/SPdeefHcxfqL1+ZGeO7Er+U8f1RVG+C5K/ZpESnp1NYZgCcBtMTzb3TmmCFAdVVdBzyEpxx8XeBmPHd3o6qHgHxe/Wj8SKjbAZiAE4WnxC94CnpF4flGDrBEVXcCiKdM8g3AUWCrqm5z2sQCZ7+lZoMA/xGRpnhKCpcBrnH2bVfVRZm8NlpVk84eyJM/JjpPb8Hzx/D/nO15gYVADWCbqm52XjM+G/HeDTQAljrHKsDfBc1SgDjn8TI8f4TBUyf/n+CpJoqnjPGfIvK7iNRzPuMKVf093XuVAv5It22JV/+C549/R+dxOaAqnoT2varudz7XRKCa06YxnnLf4KmtM1Y8xf+mqWqy13H3AdcD6WMyuZwlApNjnG+DzYEbRUTxrCqnIvKS0+SkV/M0PP//MirLe0Yq55615s+gTTSes4gGqnpKPBUlz7Q7dtEf4u/XCJ5a/ucMd4hIXS5cNvhC8Qrwmar2zeA1p/TvOi9n+iQzHwPd8PzhHpvB/uOc309n+0FEmgEt8Cya8peIfO/V/kKfqzUQD6Cq852k2xb4XETeUtX/Oe3yO+9v/IwNDZmcdD+eoZgKqnqDqpbDU13x9kxeswGo5DX2/w+vfb/gGRbBGWKqmMHri+GpLX9KRO4CKlzeRzhrEZ6hqyrO+xcUkWpOvBVFpLLTzjtRXCjeOcD9InK1s6+EiGQV5xygp9M+RESKOtun4hmvb0gGk7fOEE2IM4SVkWLAIScJ1MBz5gOeb/zNnGGiMOABr9fc7cSDE/c+Vf0Iz1DVmc8reJLTL1l8LpMLWSIwOSkKzx8qb5PxjCtnSFWP47naJN6ZpN2LZxjkzGtLOMNIPfFUnk3vCyBcRJLwnB1suJwP4BXXfjzfvGNFZBWexFBDVU/gGQr61onXe+4hw3idsfX+wGznWN/hWeYwM88Ad4nIajxDRrWdY6UA8/CUNU67wGtnc+HkGw+EOnEMdD4X6lmRKwbP8FciznCeiJTGM29yZgGgZkCyiKwA7gOGO9sb4JnrSM3ic5lcyKqPGteJSGFVPep8qxwJbFbVYW7HlR3OUMsLqhp5hd4vD54/0g+cmafIoE09oLeqPpwD79cFKKuq/82i3XBguqrOudz3NFeezRGY3OBxEemKZ0J2BZ6riEw6IlILz8Ty1AslAQBVXeFcMhqSyVlDtqjq+Gw2XWNJwH/ZGYExxgQ5myMwxpggZ4nAGGOCnCUCY4wJcpYIjDEmyFkiMMaYIPf/Aez5asEsEWvQAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# First, plot the data again.\n", "plt.errorbar(omega, Vratio, errVratio, fmt = 'ko', markersize = 5,\\\n", " linewidth = 1.8,\\\n", " markeredgecolor = 'b',\\\n", " markerfacecolor = (.49, 1, .63),\\\n", " capsize = 5)\n", "plt.xlabel('Angular Frequency (rad/s)')\n", "plt.ylabel('Voltage Ratio');\n", "\n", "# Here are some anagular frequency points that will passed to our fit model.\n", "xx = np.arange(1, 8e5, 100)\n", "\n", "# Here is the plot of the best-fit function.\n", "plt.plot(xx, LRCFunc(xx, a_fit[0], a_fit[1], a_fit[2]), 'k-')\n", "plt.axis((0, 4.6e5, 0, 0.75));" ] }, { "cell_type": "markdown", "id": "twelve-transaction", "metadata": {}, "source": [ "All of this has produced an \"unweighted\" fit to the data. To include weights, all we need to do is include another option in *curve_fit()*. Everything else is exactly the same! \n", "\n", "The new option is \"sigma\" and it is simply a list of the errors in the y-values ($\\Delta V_\\mathrm{ratio}$ in the case of this example). Note that many fitting routines require you to provide the actual weights as $1/\\sigma^2$. That is not the case here. You just have to provide the absolute $y$-uncertainties. We will also specify that the errors are absolute errors with the option \"absolute_sigma=True\".\n", "\n", "- Here's the command to execute the weighted fit." ] }, { "cell_type": "code", "execution_count": 25, "id": "thorough-johns", "metadata": {}, "outputs": [], "source": [ "a_fit, cov = curve_fit(LRCFunc, omega, Vratio, sigma = errVratio,\\\n", " p0 = start, absolute_sigma = True)" ] }, { "cell_type": "markdown", "id": "sonic-average", "metadata": {}, "source": [ "- Extract the fit parameters and their uncertainties." ] }, { "cell_type": "code", "execution_count": 26, "id": "naval-allergy", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "The best-fit parameters are: \n", " A ± ΔA = 0.7243152456809248 ± 0.003647934113877031 \n", " width ± Δwidth = 66816.88081227308 ± 623.0614968262379 1/s\n", " w0 ± Δw0 = 213875.00919865974 ± 310.60920671532733 1/s\n" ] } ], "source": [ "print('The best-fit parameters are:\\\n", " \\n A \\u00B1 \\u0394A =', a_fit[0], '\\u00B1', np.sqrt(np.diag(cov))[0],\\\n", " '\\n width \\u00B1 \\u0394width =', a_fit[1], '\\u00B1', np.sqrt(np.diag(cov))[1], '1/s'\\\n", " '\\n w0 \\u00B1 \\u0394w0 =', a_fit[2], '\\u00B1', np.sqrt(np.diag(cov))[2], '1/s')" ] }, { "cell_type": "markdown", "id": "chronic-debut", "metadata": {}, "source": [ "We can now report:\n", "\n", ">$A = 0.724\\pm 0.004$
\n", "$\\gamma=66800\\pm 600~\\mathrm{s}^{-1}$
\n", "$\\omega_0=213900\\pm 300~\\mathrm{s}^{-1}$\n", "\n", "- Plot the data and the best-fit curve." ] }, { "cell_type": "code", "execution_count": 27, "id": "christian-parade", "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEGCAYAAABo25JHAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAA5aklEQVR4nO3dd3hUZfbA8e9JASmioqg0KS4dqaEpShNCCQaCSHFdF11ZVkVQdxXRHwRQUdcGKrqogIoCshCqkEREOqEYegcFggXUFQwopJzfH3PBAUIIJMOdcj7PM09m7ty5c+Yl5Mx93/eeV1QVY4wxoSvM7QCMMca4yxKBMcaEOEsExhgT4iwRGGNMiLNEYIwxIS7C7QAu1DXXXKMVK1Z0OwxjjAkoa9eu/VFVS+X0XMAlgooVK7JmzRq3wzDGmIAiInvP9Zx1DRljTIizRGCMMSHOEoExxoQ4SwTGGBPiLBEYY0yIs0RgQlJ8fDwi4tzCEYlB5P+cn+GnnitZsiT/+te/+OWXX9wO2RifsURgQlJ8fDyqSmam0rnNEarJcP5FUaoyjKIk4JlZ3Yn//e8hXn55K1dddfWp5BAfH+9y9MYUrIC7jsCYgjRvHhxISWOjNiGSTJ5jCLVlNSWqH+LXbXu4QxNJKvIS32Q+RPFr+7F69UpKly7tdtjGFCg7IzAhLTUV2h2dQSSZAESSyU26i993/8BGbcJLDGb1b3WpGFGJQ4eiuPvuu8nKynI5amMKliUCE9Lq14ekYl3IcE6OM4hgeaGWdMo4PTnc/tsMTpzowsKFtxAREXvaOIJ1FZlA59NEICLtRWS7iOwSkUE5PP8vEVnn3DaJSJaIlPRlTMbAH4PFnTuHsy19O7VI4QmepxYp/HxiNwl6+6nk8DuFmBh+P1Wpc2ocIbr5j2RmKqpqicAEPPHVUpUiEg7sANoCacBqoJeqbjnH/p2BR1W1dW7HjYqKUqs1ZArS3Xf/hWnTjtKw4d9YvvwtIJGiTKMs5ehCIhPpzuUcZxP1iCSTDCKoF/EVLybcREyM29EbkzcislZVo3J6zpdnBI2BXaq6R1VPAJOB2Fz27wVM8mE8xpxl7969TJnyCQ89VJFlyzqgOgfVDI5k3sGrsxtQYsRTtL77T3SVOad1FXXKnMOXX/7ibvDGFBBfJoKywH6vx2nOtrOISFGgPTDtHM/3FZE1IrLm0KFDBR6oCS1ZWTBnDowYAY88kgSEM3DgwNP2CQ+HmBh45hno2ROSzxhHmEE0Gzd+dOmDN8YHfJkIJIdt5+qH6gwsU9Wfc3pSVceqapSqRpUqlWM5bWNy9ccFZOGUiJjF453X8uuQkWyb1YDCWVO54YaKZw385jaOcIA0kpIG2mCxCQq+HCNoBsSrarTz+CkAVR2Zw74JwFRV/eR8x7UxApMfc+bA0F7bWZle+1R/f1SRDTz3aY1c+/uzsjzXHKxbB/XqwY037qBmzWrEx8czdOjQSxW+MRfNrTGC1UAVEakkIoWAnsCsHIK7AmgBzPRhLMYAOV830OH3Waxbl/vrvLuKYmKgRo2qxMTEMGbMGI4fP+77wI3xIZ8lAlXNBB4GEoGtwKequllE+olIP69duwJJqnrUV7EYc1JO1w0kFetKvXoXfqyBAwdy8OBBJk+eXLBBGnOJ+axryFesa8hcjPj4eIYNGwaEUZSEU1NDZxDNAdI4Rlcgm6FDh+a5z19VqVmzJiVLlmTZsmW+DN+YfHOra8gYv3GyyJxqFoczOvNLqTeYWbUsr85uwJHMO1DNytPFYd5VS8PCIti27UaWL299VtVSG0A2gcTOCEzISUlJoWnTpkyYMIF77733oo6RlQVdo4+yY8F27iCRzyI7Ufm2G0lILEZ4eAEHbEwBsDMCY7xMnTqVyMhIYmNzu74xd6eqluIpTJea0ZC0lQeYN68AAzXmErFEYEKKqjJ16lTatWvHlVdeedHHyWn2UbtjCeedfWSMP7JEYELKqlWr2LdvH927d8/XcXKafTQrrMNFzT4yxm2WCExImTZtWr66hXK72nhv1jd07hxug8Um4NhgsQkpderUoVSpUixYsCDfx/K+2rh48V08+mg1xo59hwceeCD/gRpTwGyw2Bjg22+/ZePGjbRr165Ajud9tfGAATdSrVoVPvnkvFVSjPE7lghMyEhOTgYgOjq6wI8tIvTu3ZtFixaRlpZW4Mc3xpcsEZiQkZiYyHXXXUedOnV8cvzevXujqlZywgQcSwQmJGRnZ5OcnEzbtm0JC/PNr/2f/vQnGjdubN1DJuBYIjAhITU1lR9//NEn3UInZWVB3bpPk5raibff3k9Wls/eypgCZYnAhITExEQA2rZtW6DHPXPBm0XvluVfFOX1Bw9SImLWqfpDNp3U+DObPmpCQsuWLTl8+DCpqak+OX5OC940Kb6Z4ZOq2gL3xi/Y9FET0n799VeWLVvm026hHEtOHLWSEyYwWCIwQW/hwoVkZmYW2PUDOcmp5MScyE5WcsIEBEsEJuglJSVRtGhRbrnllgI/dm4lJ74+scdKTpiAYGMEJuhVqVKFqlWrMnfuXJ++j3fJiW+//Yy33+7Mxo3rqV27tk/f15i8sDECE7L27NnDrl27fDo+cJJ3yYmhQxsSFgaffvqpz9/XmPyyRGCCWlJSEuCbshK5ue6662jRogVTp04l0M66TejxaSIQkfYisl1EdonIoHPs01JE1onIZhFZ5Mt4TOhJSkrihhtuoGrVqpf8ve+66y62bdvGpk2bLvl7G3MhfJYIRCQceAvoANQEeolIzTP2uRIYA9yhqrWA/K0WYoyXjIwMFixYQHR0NCJyyd8/Li6OsLAw6x4yfs+XZwSNgV2qukdVTwCTgTNXA+kNTFfVfQCqetCH8ZgQk5KSwpEjR3w6bTQ31157La1ateLTTz+17iHj13yZCMoC+70epznbvFUFrhKRL0VkrYj8JacDiUhfEVkjImsOHTrko3BNsElKSiIsLIw2bdq4FsNdd93Fjh072LBhg2sxGHM+vkwEOZ2Ln/m1KAJoCHQCooH/E5GzOnNVdayqRqlqVKlSpQo+UhOUEhMTadKkCVdddZVrMXTt2pXw8HCmTp3qWgzGnI8vE0EaUN7rcTng2xz2ma+qR1X1R2AxUNeHMZkQ8fPPP7N69WrXuoVOKlWqlHUPGb/ny0SwGqgiIpVEpBDQE5h1xj4zgVtFJEJEigJNgK0+jMmEiM8//xxVveTTRnNy1113sXPnTtavX+92KMbkyGeJQFUzgYeBRDx/3D9V1c0i0k9E+jn7bAXmAxuAVcB7qmpz7Uy+JSYmcsUVV9CoUSO3QznVPWSzh4y/shITJmicLPHw1VfK6NF9aNHiGNOm+ccf3+joaHbv3s3OnTtdmcpqjJWYMEHrzIVhHu+8lvShL3D1T/2ZP/3PfrMwTPfu3dm9e7fP1kMwJj8sEZiAFh8fj6oye3YW1YtXYxNNeYnBbKIpVYtWZfbsLFTV9URg3UPGn1kiMEEhp4Vh2v82028Whrn66qu5/fbbmTJlis0eMn7HEoEJCjktDJNYrKtfLQzTu3dvvvnmG1auXOl2KMacxhKBCWi5LQyzPX2bXy0M06VLFy677DImTZrkdijGnMZmDZmgkZUF3buPY+bMvUyePIi4uCKEh7sd1em6d+/O4sWLOXDgABEREW6HY0KIzRoyISE8HHbvHkWLFkvo3t3/kgBAr169OHjwR0aO3MiIETBnjieBGeMmSwQmaHz33Xds2LDBL64mPtPJLqxu3bpTlAQmDsnm1yEjebzzWkpEzPKbaa4mNFkiMEEjOTkZwPX6QjnxnuZaMaLCadNcqxWv7jfTXE1oskRggkZiYiLXXnstdev6b93C1FTonDX/tGmu0UcT/GaaqwlNlghMUMjOziY5OZm2bdsSFua/v9b160Oyn09zNaHHf//HGHMB1q1bx6FDh/xyfAACa5qrCT2WCExQSExMBKBt27YuR5Kzk2MEqlkcybyDvv8uwr85Ssyj/+NI5h2o2hiBcY8lAhMUkpKSqFu3Ltdff73boZxXeDg8/nh1Klb8mC1b/u2X01xNaLFEYAJeeno6y5Yt89tuoZyICHfffTfJyckcOHDA7XBMiLNEYALewoULycjI8Mtpo7m59957yc7OZuLEiW6HYkKcJQIT8BITEylatCjNmzd3O5QLUqVKFZo3b86ECROsIqlxlSUCE/Dmz59P69atKVy4sNuhXLC//vWvbNu2jVWrVrkdiglhlghMQNu1axe7d++mffv2bodyUbp3706RIkWYMGGC26GYEGaJwAS0k9NGAzURlChRgm7dujFp0iR+//13t8MxIcqniUBE2ovIdhHZJSKDcni+pYgcFpF1zm2IL+MxwWf+/PnceOON3HjjjW6HctH++te/cvjwYWbOnOl2KCZE+SwRiEg48BbQAagJ9BKRmjnsukRV6zm34b6KxwSf48eP88UXXwTs2cBJrVq1onz58tY9ZFzjyzOCxsAuVd2jqieAyUCsD9/PhJilS5dy7NixgE8EYWFh3HvvvSQlJZGWluZ2OCYE+TIRlAX2ez1Oc7adqZmIrBeReSJSK6cDiUhfEVkjImsOHTrki1hNAJo/fz6FChWiZcuWboeSb/fddx+qyvvvv+92KCYE+TIRSA7bzpws/RVQQVXrAm8AM3I6kKqOVdUoVY0qVapUwUZpAlZiYiK33norxYsXdzuUfKtUqRLt2rXj3XffJTMz0+1wTIjxZSJIA8p7PS4HfOu9g6oeUdV05/5nQKSIXOPDmEyQOHDgABs3bgz4biFv/fr148CBA8ydO9ftUEyI8WUiWA1UEZFKIlII6AnM8t5BRK4XEXHuN3bi+cmHMZkgEejTRnMSExNDmTJleOedd9wOxYQYnyUCVc0EHgYSga3Ap6q6WUT6iUg/Z7c7gU0ish4YDfRUu9be5MH8+fMpW7YstWrlOKwUkCIiIvjb3/5GYmIiX3/9tdvhmBAiefm7KyJ3ALc5Dxep6myfRpWLqKgoXbNmjVtvb/xARkYGpUqV4s477+S9995zO5wCtX//fipWrMiTTz7J888/73Y4JoiIyFpVjcrpufOeEYjISGAAsMW5PeJsM8YVS5Ys4fDhw8TExLgdSoErX748MTExvP/++5w4ccLtcEyIyEvXUCegraqOU9VxQHtnmzGumD17NoULF/bb1cjyq1+/fhw8eJBp06a5HYoJEXkdI7jS6/4VPojDmDxRVWbPnk3r1q0pVqyY2+H4RHR0NFWqVGHUqFFuh2JCRF4SwUggVUQmiMgHwFrAOi+NK7Zv387u3bvp3Lmz26H4TFhYGAMGDCAlJYUVK1a4HY4JAedNBKo6CWgKTHduzVR1sq8DMyYns2d75ikE4/iAt3vvvZcrrriC119/3e1QTAg4ZyIQkerOzwZAaTwXiO0HyjjbjLlksrJgzhx4662rqFTpYcqUKX/+FwWw4sWL88ADDzBt2jT27dvndjgmyOV2RvCY8/OVHG4v+zguY07JyoKu0UcZ0nMbd+09RLG0vnSNPkpWltuR+Vb//v0BePPNN12OxAS7cyYCVe3r3O2gqq28b0DHSxOeCWXx8fGICBERMWxfsI2UozfxEoP5KqMB2xZsJyIiBhEhPj7e7VB94oYbbiAuLo53332X9PR0t8MxQSwvg8XL87jNmAIVHx+PqjJ8+Bzi5HMi8RRjiySTbpLMiBFzUNWgTQQAjz76KL/88otVJTU+ldsYwfUi0hAoIiL1RaSBc2sJFL1UARpTvz4kFe1CBhEAZBBBYrGu1KvnblyXQrNmzbj11lt5+eWX7QIz4zO5nRFE4xkLKAe8yh/jA48Bg30fmgl1J7uGOncOZ9vR7dQihSd4nlqksD19G507hwd119BJgwcPJi0tjYkTJ7odiglS5601JCLdVNVvLnG0WkOhqU+fvzFlyhGeeOIToqIi6NABwsPdjurSUFUaNmxIeno6W7duJTxUPrgpULnVGoo434tVdZqIdAJqAZd5bbf1hc0lkZmZyezZM+jaNZr4+PP+ygYdEWHw4MF0796dadOmcdddd7kdkgkyeSk69w7QA+iPZ9Wx7kAFH8dlzCmLFi3ip59+olu3bm6H4pq4uDiqV6/O888/j1VqNwUtL7OGblbVvwD/U9VhQDNOX3nMGJ+aPn06RYsWDapFaC5UWFgYgwYNYv369baCmSlweUkEvzk/j4lIGSADqOS7kIz5Q3Z2NgkJCXTo0IGiRUN7slrv3r2pXLkyQ4YMITs72+1wTBDJSyKYIyJXAv/Gs9j8N4DVGjKXxIoVK/juu++Ii4tzOxTXRUZGEh8fT2pqKtOmTTtVdmPECM/PYL/S2vhOnlYoO7WzSGE8A8aZqnrUZ1HlwmYNhZbHHnuMt956i0OHDlGiRAm3w3FdVlYWderUISsLqpRN4dtVB2h3dAZJxbpQtkk5EhKLhcxsKnNhLnqFMhEpKyJRzuLz4FmL4ElgZwHHaMxZsrOzmTZtGm3btrUk4BgxYgRbtmxh+/ZK7PhiOyvTazNSB7EyvXZIlN0wvpHblcUDgXXAG8BKEbkXzyL0RYCGlyI4E9qWL1/Ovn376NWrl9uh+I34+Hiys7MpXboTsSSFZNkNU/ByOyPoC1RT1WZAF+BdoJOqPqqq3+Xl4CLSXkS2i8guERmUy36NRCRLRO68kOBNcPvkk08oUqQIsbGxbofiV0SEvn0bMZN2IVl2wxS83BLB76r6M4Cq7gN2qOrKvB5YRMKBt4AOQE2gl4jUPMd+LwKJFxK4CW4ZGRlMnTqVO+64g+LFi7sdjt84WXZj2LAmpHEgZMtumIKV22Wa5URktNfja70fq+oj5zl2Y2CXqu4BEJHJQCyw5Yz9+gPTgEZ5jtoEvc8//5wff/zRuoXOEB8ff+qP/Jo1qTRqNITltzzIq4Ma0KFDA8LDbeqQuXC5JYJ/nfF47QUeuyyeFc1OSgOaeO8gImWBrkBrckkEItIXT1cVN9xwwwWGYQLRpEmTuPLKK0P6IrLziYqqT58+pZg4MZbq1bcQHv4nt0MyAeqciUBVP8jnsSWnw57x+HXgSVXNEslp91OxjAXGgmf6aD7jMn7u2LFjJCQk0KNHDwoXLux2OH7tueee49NPP+XJJ59k2jS/qQ1pAkxeLii7WGmcXoqiHPDtGftEAZNF5BvgTmCMiHTxYUwmAMyZM4f09HR69+7tdih+r3Tp0gwaNIjp06ezaNEit8MxAeqCLii7oAOLRAA7gDbAAWA10FtVN59j/wnAHFX9b27HtQvKgl9sbCyrV69m//79VnI5D44dO0aNGjW4/PLLSU1NJTIy0u2QjB+66AvK8kNVM4GH8cwG2gp8qqqbRaSfiPTz1fuawPb9998zd+5c7rnnHksCeVS0aFHeeOMNNm/ezGuvveZ2OCYA5WVhmqrA28B1qlpbROoAd6jqs5ciwDPZGUFwe/nll/nXv/7F1q1bqV69utvhBJQuXbqQlJTEli1bqFixotvhGD+T3zOCd4Gn8FQdRVU3AD0LLjxjPFSV8ePH07RpU0sCF2H06NGICI88cr6Z3cacLi+JoKiqrjpjW6YvgjGhbfXq1WzZsoU+ffq4HUpAuuGGGxg2bBizZ89mxowZbodjAkheEsGPInIjztRPpwxEnkpMGHMhxo8fT5EiRejRo4fboQSsAQMGUKdOHR588EF+/vlnt8MxASIvieAh4D9AdRE5AAwE/uHLoEzo+e2335g0aRLdunXjiiuucDucgBUZGcn48eM5dOgQAwYMcDscEyDOmwhUdY+q3g6UAqqranNV/cbnkZmQkpCQwOHDh61bqAA0aNCAp59+mokTJ1oXkcmTvMwaeiyHzYeBtaq6zhdB5cZmDQWnFi1akJaWxs6dOwkL8+V1jqHhxIkTNGnShG+//ZbNmzdzzTXXuB2ScVl+Zw1FAf3w1A4qi6fmT0vgXRF5oqCCNKFr06ZNLF68mH79+lkSKCCFChViwoQJ/O9//+PBBx/EVxeOmuCQl/91VwMNVPVxVX0cT2IoBdwG/NWHsZkQ8fbbb1O4cGHrFipgdevWZfjw4UydOpVx48a5HY7xY3lJBDcAJ7weZwAVVPU34LhPojIh49dff+Wjjz7irrvusu4LH3jiiSdo06YN/fv3Z+vWrW6HY/xUXhLBJ3iWqhwqIkOBZcAkESnG2WsLGHNBPv74Y3799VcefPBBt0MJSmFhYXz00UcUL16cnj178vvvv7sdkvFDeSo6JyJRwC14SksvVVXXRmttsDh4qCr16tUjLCyMr776itxKkZv8mTdvHh07duTBBx/krbfecjsc44J8F51z/vBPAqYDB0XEVocxFy0rC+bMgQce2MeGDeXp1+8hSwI+1qFDBx5//HHGjBnDxIkT3Q7H+JnzJgIRuUNEdgJfA4ucn/N8HZgJLifX2hUJp0TELB7vvJaS739CVYbxWL9rEbG1dn1t5MiRtGjRggceeIDU1FS3wzF+JC/XEazHs5Tk56paX0RaAb1Ute+lCPBM1jUU2ObMgaG9trMyvTaRZJJBBE2Kb2b4pKrExLgdXfA7ePAgDRs2JDw8nDVr1tgAfQjJb9dQhqr+BISJSJiqLgTqFWSAJnSkpkK7ozOIdOoWRpJJ9NEE1q1zN65Qce2115KQkMD3339Pjx49yMy0+pEmb4ngFxEpDiwGPhaRUVj1UXOBTnYNDRkSw3S9nQxnuewMIpimbfm//4uxrqFLJCoqinfeeYcvvviCRx55hMxMZc4cGDHCc8aWleV2hOZSy0vXUDHgNzxJ427gCmCiqrpS2tC6hgJbVhbUr7qZ3/ccJ06SSSrWlXJNypKQWAxbkOzSiY+PZ9iwYUAYRUmgHGWJJYmZtCONAxyjK5DN0KFDLTkHidy6hvKSCF5U1SfPt+1SsUQQ2I4fP06FCpUpU+Z+4uKGU68edOiAJQEXZGdn06LFvzm49HY20dTGbIJcfscI2uawrUP+QjKh6oMPPuCHH77lxRdv5ZlnICbGkoBbwsLCaN36MbqQbGM2Ie6ciUBE/iEiG4FqIrLB6/Y1sOHShWiCRUZGBiNHjqRJkybcfvvtbocT0k6O2Qwf3pUZtLUxmxB3zq4hEbkCuAoYCQzyeurXvI4PiEh7YBQQDrynqi+c8XwsMALIxjMAPVBVl+Z2TOsaClzjx4/nvvvuY86cOXTq1MntcAyeMZuu0UfZv2I/bY8lMJP2lG9WjsQlpexMLchc1BiBiJTM7aDnSwYiEg7swNO1lAasxnP9wRavfYoDR1VVRaQO8Kmq5rpquSWCwJSZmUmNGjUoUaIEa9assSuJ/UhWFsybB8nJh5gwYSAlSixj6dJFVKhQwe3QTAHKLRFE5PK6tTjrFOOpMeRNgcrned/GwC5V3eMEMRmIxatQnaqme+1fzOv9TJCZMmUKu3btYvr06ZYE/Ex4uGesJiamFPfd9wQtW7akTZs2LFq0iLJly7odnrkEzjlGoKqVVLWyc6t0xu18SQA8i9js93qc5mw7jYh0FZFtwFzgvpwOJCJ9RWSNiKw5dOhQHt7a+JOsrCyeffZZateuTWxsrNvhmFzUrVuXefPm8cMPP9CiRQv27t3rdkjmEshT0Tmn3tDLzi2vk8py+tp31jd+VU1wuoO64BkvOPtFqmNVNUpVo0qVKpXHtzf+4uOPP2bbtm0MGTLEViALAE2bNiU5OZkff/yR2267jd27d7sdkvGxvBSdewEYgKdLZwswQERG5uHYaUB5r8flgG/PtbOqLgZuFBErfhJEjh8/zpAhQ2jQoAHdunVzOxyTR02bNuWLL77g6NGj3HbbbWzbts3tkIwP5eXrWUegraqOU9VxQHsgL1M+VgNVRKSSiBQCegKzvHcQkT+J02EsIg2AQsBPF/IBjH/7z3/+w969e3nhhRfsbCDANGjQgC+//JLMzExatGiBTdIIXnn9n3ml1/0r8vICVc0EHgYSga14ZgRtFpF+ItLP2a0bsElE1gFvAT3UVtkOGr/++ivPPvssrVu3tusGAlTt2rVZvHgxRYsWpWXLlsybZxXog1Fus4ZOGgmkishCPP3+twFP5eXgqvoZ8NkZ297xuv8i8GKeozUB5dVXX+XQoUOMHDnSZgoFsGrVqrFixQo6duxI586d+c9//sP999/vdlimAOV2ZfGbInKzqk4CmuJZnWw60ExVJ1+qAE1g+v7773n55ZeJi4ujcePGbodj8un6669n0aJF3H777fztb39j6NChZGdnux2WKSC5dQ3tBF4RkW+AgcA+VZ2pqt9fisBMYBs8eDDHjx/nhRdeOP/OJiBcfvnlzJ49mz59+jB8+HC6d+9Oenr6+V9o/F5u1xGMUtVmQAvgZ2C8iGwVkSEiUvWSRWgCzurVqxk/fjyPPvooVapUcTscU4AiIyN5//33efXVV5kxYwbNmjVjz549bodl8um8ZahP21mkPjAOqKOqrlQisRIT/k1VueWWW9izZw87duygRIkSbodkfCQ5OZkePXogIkyZMsUmBPi5fJWhFpFIEeksIh/jWbR+B57ZPsac5ZNPPmHFihWMHDnSkkCQa9u2LatXr6Z06dJER0czbNgwsmx5s4CUW9G5tkAvPNcMrAImAzNU9eilC+9sdkbgv44cOUKNGjUoU6YMKSkpdt1AiEhPT+ehhx7iww8/pGXLlnz88ceUKVPG7bDMGS72jGAwsAKooaqdVfVjt5OA8W/PPPMM3333HW+++aYlgRBSvHhxPvjgAyZMmMCqVauoV68e8+fPBzyVTW09ZP93QWME/sDOCPzTqlWraNq0KQ899BBvvPGG2+EYl2zZsoVWrVpx8OBBbD1k/5KvNYv9jSUC/5ORkUFUVBQ//fQTW7ZssbGBEHfs2DEGDx7MqFG7qCbD2ahNbD1kP5DfNYuNydVrr73Ghg0beOONNywJGIoWLcrrr79Onz5vcIcm2nrIAcASgcmXHTt2EB8fT2xsLF27dnU7HOMHTq6HPH58f2bSztZDDgDWNWQuWmZmJs2bN2fHjh1s2rTJZoqY05xcDzkt5QDtjiYwg3bs1/106f1fXnvtZa699lq3Qwwp1jVkfOKll14iJSWFt956y5KAOUt4OCQkFmP4pKoUH/4kz02tycDBa5g6dTLVqlVj7Nixdt2Bn7AzAnNR1q9fT6NGjejatSuTJ0+26qImz7Zu3co//vEPFi1aRL169Xjttddo2bKl22EFPTsjMAXq+PHj3HPPPVx99dWMGTPGkoC5IDVq1GDhwoV88skn/PTTT7Rq1Yq4uDh27drldmghyxKBuWBPPPEEGzdu5L333uPqq692OxwTgESEXr16sX37dp599lmSkpKoWbMmjz/+OD/9ZIsUXmqWCMwFSUhIYPTo0QwYMIBOnfKyYqkx51akSBGefvppdu7cyT333MNrr71GpUqVGDp0KIcPH3Y7vJBhicDk2TfffMN9991HVFQUL730ktvhmCBSunRp3n//fTZu3Ei7du0YPnw4lSpV4oUXXuDoUats42uWCEyenDhxgp49e5Kdnc2UKVMoVKiQ2yGZIFSrVi3++9//snbtWpo1a8ZTTz1F5cqVefHFF+0MwYcsEZg8eeKJJ0hJSeH999+ncuXKbodjglyDBg2YO3cuy5Yto27dugwaNIgKFSowePBgp46RKUg+TQQi0l5EtovILhEZlMPzd4vIBue2XETq+jIec3E++OADRo0axcCBA7nzzjvdDseEkJtvvpmkpCRWr15N27ZteeGFF6hQoQIPP/wwX3/9tdvhBQ9V9ckNCAd2A5WBQsB6oOYZ+9wMXOXc7wCknO+4DRs2VHPprFy5UgsXLqxt2rTRjIwMt8MxIW7btm16//33a2RkpIaFhWmXLl30iy++0OzsbLdD83vAGj3H31VfnhE0Bnap6h5VPYFnYZvYM5LQclX9n/NwJVDOh/GYC/Tdd98RFxdHmTJlmDJlChEREW6HZEJctWrVeO+99/j6668ZNGgQS5YsoXXr1tSpU4d3332XY8eOuR1iQPJlIigL7Pd6nOZsO5f78SyFeRYR6Ssia0RkzaFDhwowRHMuv/32G3FxcRw+fJiZM2fa9QLGr5QtW5bnnnuO/fv3M27cOCIiIujbty/lypXjscceY8uWLW6HGFB8mQhyutw0x3oWItIKTyJ4MqfnVXWsqkapalSpUqUKMESTk6ysLO6++y+sXHkNnTuvYu/em2xlKeOXihQpQp8+ffjqq69YsmQJbdu25c0336RWrVrcfPPNjBs3jvT09NNeY6umnc2XiSANKO/1uBzw7Zk7iUgd4D0gVlXtkkI/8Oij/yQx4R5uKjSSilNmM7TXdrpGH7X/MMZviQjNmzdnypQpHDhwgFdeeYVffvmF+++/n9KlS9O3b19WrlxJZqbSNfooQ3tt59jQF+13+6RzDR7k9wZEAHuASvwxWFzrjH1uAHYBN+f1uDZY7DtDhw5VQKGTVmWNniBCFfQEEVqFtQqdFNChQ4e6Haox55Wdna3Lli3TPn36aGRkZMj/bpPLYLHPRv9UNVNEHgYS8cwgGqeqm0Wkn/P8O8AQ4GpgjFO4LFPPUR3P+F6NGjUAqFnzz3Te+jmR+sfKUt0kmWLD5/DMM25GaEzeiQg333wzN998M6+//jrTp09nxAiI3ZN02qppcZJM8VD/3T5XhvDXm50R+MasWbM0LCws5L81meCU17PdgQMHuh2qz5DLGYGtR2BITk4mJiaGevXqkZycTLFiJU6tLBV9NIHEYl0p16QsCYnFCA93O1pjLp73qmnRRxOYVySWzOsykGL3sGnTegCioqKIi4ujW7duVK1a1eWIC05u6xFYIghxS5YsITo6mipVqrBw4UJKliwJeP7DzJsH69ZBvXrQoQOWBExQONfv9s6dO0lISGD69OmkpKQAntpHsbGxdOzYkaZNmxIewP8JLBGYHC1fvpz27dtTpkwZFi9ebGvIGuNIS0tjxowZTJs2jSVLlpCVlUXJkiWJjo6mY8eOtG/fnmuuucbtMC+IJQJzli+//JKYmBhKly7NwoULKVfOLuo2Jie//PILycnJzJ07l3nz5nHw4EFEhCZNmtCxY0fatm1LVFSU3195b4nAnCYpKYnY2FgqV67M559/TunSpd0OyZiAkJ2dzVdffcVnn33G3LlzWb16NapKiRIlaNGiBW3atKFNmzbUqlXL75ZwtURgTpk9ezZ33nknNWrUIDk5GbtS25iLd+jQIRYuXMiCBQtYsGABu3fvBuC6666jdevWtGnThlatWlGpUiXXE4MlAgPAhx9+yP3330/9+vVJTEzkqquucjskY4LK3r17TyWFBQsW8MMPPwBQpkwZmjdvTvPmzbn11lu56aabLvnAsyWCEKeqjBw5kqeffpo2bdowffp0SpQo4XZYxgQ1VWXLli0sXryYpUuXsmTJEvbv99ThvPzyy7n55pu59dZbad68OY0aNaJo0aI+jccSQQjLysqif//+vP3229x9992MGzfOlpk0xiX79u07lRSWLl3Kpk2bAAgPD+emm26icePGNGnShMaNG1OjRo0CPWuwRBCijh07Ru/evZk5cyZPPPEEI0eOJCzMVic1xl/8/PPPLF++nJUrV7Jq1SpWrVp1am3mYsWKERUVRZMmTWjYsDFHj97G/v3X0KCBXNR1PZYIQtC+ffuIjY1l/fr1jBo1iv79+7sdkjHmPLKzs9m1axcpKSmnEkNq6gYiM6ZQjrLEksSc8A5QVnj6+a00bFiPKlWq5OnMIbdE4N8TX81FWbp0KXFxcRw/fpw5c+bQsWNHt0MyxuRBWFgYVatWpWrVquzevZtVq1YBnahEWTbRlEgyeS5rCLX2pfDnP08EehEZGUnDhg2pX78+9evXp169etSuXZsiRYrk+X3tjCDIvPfeezz44INUrFiRWbNmUb16dbdDMsbkw4gRcGzoi4zUQae2DZIXSH/wL0RFJbJu3TpSU1NZt24dR44cATwJpUqVKtx0002nbnFxcXZGEOyOHz/OY489xpgxY4iOjmbSpEk2PdSYABYfH8+wYcOATlRlGMOJIJJMMohgurZl51sPAHMZOnQoixYtIjs7m2+++YbU1FTWr1/Pxo0bSU1NZdq0aZzvC7+dEQSBPXv2cNddd7F27Vr++c9/MnLkSL+/3N0YkzdnVky90GrAR48eZfPmzTRp0sQGi4NRVhYMGbKCV175goiIjXz4YS/i4mLdDssYU8AKohqwzRoKQr/9doL61Taj+7PpQjJJRbtQvll5WzPAGJOj3BKBTSoPMPHx8YgIRYvGofuz2URTXuQpVh27iW0LthMREYOIEB8f73aoxpgAYR3JASQ7O5srr7ySwoULEx7enK6/2brCxpj8szOCAJGWlka7du149NFHOX78OMeOLSVBbyfDyeUZRDBN2/J//2dnBMaYC+PTRCAi7UVku4jsEpFBOTxfXURWiMhxEfmnL2MJVKrKBx98wE033cTKlSsZO3Ys2dnZZGbOoVqb6jQpvpmn5AWaFN9M9TbVyMycg6paIjDG5JnPuoZEJBx4C2gLpAGrRWSWqm7x2u1n4BGgi6/iCGRff/01f//730lOTqZ58+aMHz+eP/3pT4BnxkBCYjHmzavKunVPMryerStsjLk4vhwjaAzsUtU9ACIyGYgFTiUCVT0IHBSRTj6MI+BkZWUxevRonnnmGcLDwxkzZgx///vfzyoYFx4OMTGemzHGXCxfJoKywH6vx2lAk4s5kIj0BfoC3HDDDfmPzI9t3LiRBx54gJSUFGJiYhgzZgzly5d3OyxjTBDz5RhBTuuyXdRFC6o6VlWjVDUqWJdWPHz4MAMHDqR+/frs3r2bSZMmMWvWLEsCxhif82UiSAO8/4qVA7714fsFpOzsbD744AOqVq3K6NGj6du3L9u3b6dnz56ur3FqjAkNvuwaWg1UEZFKwAGgJ9Dbh+8XcNatW8dDDz3E8uXLadq0KfPmzaNBgwZuh2WMCTE+OyNQ1UzgYSAR2Ap8qqqbRaSfiPQDEJHrRSQNeAx4RkTSRCToF9P9/vvv6devHw0bNmTnzp2MGzeOZcuWWRIwxrjCp1cWq+pnwGdnbHvH6/73eLqMQsLRo0d59dVXefHFFzl+/Dj9+/cnPj6eK6+80u3QjDEhzEpM+FhWFsyZk8WECetZvHgUP/88kW7dujJy5EiqVKnidnjGGGOJwJcyM5UWjfbz84af6JydyLawgdRqOJIpU8rYhV/GGL9htYZ8ZPHixdStO4gf1x1iQ3ZjXmIwG7Ibk749nXnz3I7OGGP+YImggK1cuZIbb7yRFi1asGVLIWJJIpI/KoTenj6Nzp2fscJwxhi/YYmggKSmptK5c2eaNWvG999/f3IrM2l3WoXQGUQD69wK0xhjzmKJIJ/WrVvHnXfeSYMGDVi6dCnPP/88P/zwA6pqFUKNMQHBlqq8SCtWrOC5555j7ty5lChRggEDBvDYY4+dNRW0INYaNcaY/MptqUqbNXQBVJWFCxfy7LPPsnDhQq6++mpGjBjBww8/fM5rAaxCqDHG31kiyANVZe7cuTz33HOsXLmS0qVL88orr9C3b1+KFy/udnjGGJMvlghykZGRwdSpU/n3v//NunXrqFChAmPGjKFPnz5cdtllbodnjDEFwhJBDo4cOcK7777LqFGj2L9/P9WrV2fChAn07t2byMhIt8MzxpgCZYnAy/79+xk9ejRjx47lyJEjtGzZkjFjxtCxY8ezVgczxphgYYkAzxTQV155hcmTJ6OqdO/enccff5yoqBwH2I0xJqiEbCJQVRITE3n55ZdZsGABxYsX5+GHH2bAgAFUrFjR7fCMMeaSCfpEcHIef2oq1K8PzZv/yscff8gbb7zB9u3bKVOmDC+++CJ9+/a1ctDGmJAU1IkgKwu6Rh/lQEoabY/OYFBEB/Zlp/Fr1iM0atSQDz/8kB49elCoUCG3QzXGGNcEdSKYOzebfcv2svr3ukSSyYiMZ6gbsZb7X97C449Xczs8Y4zxC0E5FaZ58+aICLGxQ2n3+8zTqn/GZM7ln//8CBGhZcuW7gZqjDF+IKgSwY4dO3jkkUfYsGEDAFWrHiOx8B2nVf/8vHg3Zs9+FlXlyy+/dDFaY4zxDwGfCLKzs5k3bx4dOnSgWrVqvPPOO8TGxpKSksKWLa9QoXnF06p/lmtSlg4d3I7aGGP8h0+rj4pIe2AUEA68p6ovnPG8OM93BI4Bf1XVr3I75snqo0eOHGHChAm8+eab7Ny5k+uvv55//OMf9O3bl+uvv/7U/lb90xhjcq8+6rNEICLhwA6gLZAGrAZ6qeoWr306Av3xJIImwChVbZLbcWvXrq2tWrViwoQJpKen07RpUx555BG6detms3+MMeYc3CpD3RjYpap7nCAmA7HAFq99YoEP1ZONVorIlSJSWlW/O9dBN2/ezM6dO+nRowf9+/enUaNGPvwIxhgT/HyZCMoC+70ep+H51n++fcoCpyUCEekL9HUeHj9x4sSmjz76iI8++qhgIw4e1wA/uh2En7M2Oj9ro9wFWvtUONcTvkwEksO2M/uh8rIPqjoWGAsgImvOdXpjPKyNzs/a6PysjXIXTO3jy1lDaUB5r8flgG8vYh9jjDE+5MtEsBqoIiKVRKQQ0BOYdcY+s4C/iEdT4HBu4wPGGGMKns+6hlQ1U0QeBhLxTB8dp6qbRaSf8/w7wGd4ZgztwjN9tE8eDj3WRyEHE2uj87M2Oj9ro9wFTfv49DoCY4wx/i/gryw2xhiTP5YIjDEmxAVUIhCR9iKyXUR2icggt+MpaCIyTkQOisgmr20lRSRZRHY6P6/yeu4ppy22i0i01/aGIrLReW60U8oDESksIlOc7SkiUtHrNfc677FTRO69RB/5gohIeRFZKCJbRWSziAxwtlsbOUTkMhFZJSLrnTYa5my3NjqDiISLSKqIzHEeh24bqWpA3PAMOO8GKgOFgPVATbfjKuDPeBvQANjkte0lYJBzfxDwonO/ptMGhYFKTtuEO8+tAprhuU5jHtDB2f4g8I5zvycwxblfEtjj/LzKuX+V2+2RQ/uUBho49y/HU8KkprXRaW0kQHHnfiSQAjS1NsqxrR4DPgHmOI9Dto1c/8e4gH+0ZkCi1+OngKfcjssHn7MipyeC7UBp535pYHtOnx/P7Kxmzj7bvLb3Av7jvY9zPwLPVZHivY/z3H/w1IVyvT3O01Yz8dSysjbKuX2KAl/huaLf2uj0tikHLABa80ciCNk2CqSuoXOVowh216lzbYXz81pn+7nao6xz/8ztp71GVTOBw8DVuRzLbzmn2vXxfOO1NvLidHmsAw4CyapqbXS214EngGyvbSHbRoGUCPJUjiKEnKs9cmuni3mN3xGR4sA0YKCqHslt1xy2BX0bqWqWqtbD8623sYjUzmX3kGsjEYkBDqrq2ry+JIdtQdVGgZQIQrUcxQ8iUhrA+XnQ2X6u9khz7p+5/bTXiEgEcAXwcy7H8jsiEoknCXysqtOdzdZGOVDVX4AvgfZYG3m7BbhDRL4BJgOtRWQiodxGbvdNXUCfXgSegZVK/DFYXMvtuHzwOSty+hjBvzl9AOsl534tTh/A2sMfA1ir8QwQnhzA6uhsf4jTB7A+de6XBL7GM3h1lXO/pNttkUPbCPAh8PoZ262N/miLUsCVzv0iwBIgxtronO3Vkj/GCEK2jVz/h7jAf7SOeGaK7AaedjseH3y+SXhKcGfg+eZwP55+xQXATudnSa/9n3baYjvObAVnexSwyXnuTf64gvwyYCqekh6rgMper7nP2b4L6ON2W5yjfZrjOY3eAKxzbh2tjU5rozpAqtNGm4AhznZro5zbqyV/JIKQbSMrMWGMMSEukMYIjDHG+IAlAmOMCXGWCIwxJsRZIjDGmBBnicAYY0KcJQJT4ESkq4ioiFT34Xuk5/P1XzqVJNc5tzsLKja3iUgREVkkIuH5OMZfReRNr8elRSQpl/0/967WaQKLJQLjC72ApXgupHGdeOT0u363qtZzbv894zUX/UfUD9wHTFfVLO+N+fxM7fEUUjuXj/BU3DQByBKBKVBOHaBb8FwM19Nre0vnW/h/RWSbiHzsVbu9o7NtqVPT/WR9+HgR+afXMTZ513U/+X4iskBEvnLqwsc62yuKZ92CMXgqcHpf1n+u2L8RkSEishToLiLtRGSFc+ypzmc7uS7GBcUrIn8WzzoB60TkPyf/KItIuog8J571A1aKyHXO9utEJMHZvl5EbhaREeKsweDs85yIPJLDR7kbT2XWk+2+UEQ+ATY622aIyFrxrFfQ1+t4fURkh4gscv4NvbUH5jlnBoudz7FJRG51np+F5wuACUCWCExB6wLMV9UdwM8i0sDrufrAQDz13SsDt4jIZXhK8XZQ1eZ4SiRciN+BrqraAGgFvHIywQDVgA9Vtb6q7s3htR97dQ1dffJ4ThyfA88AtzvHXgM85sT7LtAZuBW4/nwBikgNoAdwi3qKwWXh+WMNUAxYqap1gcXAA8720cAiZ3sDYDPwPnCvc8wwPIn24zPeqxCeq1i/8drcGM+V+DWdx/epakM8V8U+IiJXO7V1huFJAG3x/BudPGY4UE1VtwC98ZSDrwfUxXN1N6r6P6CwVzuaABLhdgAm6PTCU+IXPAW9euH5Rg6wSlXTAMRTJrkikA7sUdWvnX0mAae+peaBAM+LyG14SgqXBa5znturqitzee3dqrrm1IE8+WOK87Apnj+Gy5zthYAVQHXga1Xd6bxmYh7ibQM0BFY7xyrCHwXNTgBznPtr8fwRBk+d/L+Ap5oonjLGh0XkJxGp73zGVFX96Yz3ugb45Yxtq7zaFzx//Ls698sDVfAktC9V9ZDzuaYAVZ19muAp9w2e2jrjxFP8b4aqrvM67kGgDHBmTMbPWSIwBcb5NtgaqC0iimdVORWRJ5xdjnvtnoXn9y+nsrwnZXL6WetlOexzN56ziIaqmiGeipIn9zt6wR/ij9cInlr+p3V3iEg9zl02+FzxCvCBqj6Vw2sy9I86LyfbJDfvAX/F84d7XA7P/8bZ7XSqHUSkJXA7nkVTjonIl177n+tzdQDmA6jqYifpdgI+EpF/q+qHzn6XOe9vAox1DZmCdCeerpgKqlpRVcvjqa7YPJfXbAMqe/X99/B67hs83SI4XUyVcnj9FXhqy2eISCugQv4+wikr8XRd/cl5/6IiUtWJt5KI3Ojs550ozhXvAuBOEbnWea6kiJwvzgXAP5z9w0WkhLM9AU9/fSNyGLx1umjCnS6snFwB/M9JAtXxnPmA5xt/S6ebKBLo7vWaNk48OHEfVNV38XRVnfy8gic5fXOez2X8kCUCU5B64flD5W0ann7lHKnqb3hmm8x3Bml/wNMNcvK1JZ1upH/gqTx7po+BKBFZg+fsYFt+PoBXXIfwfPOeJCIb8CSG6qr6O56uoLlOvN5jDznG6/StPwMkOcdKxrPMYW4GAK1EZCOeLqNazrFOAAvxlDXOOsdrkzh38p0PRDhxjHA+F+pZkSseT/fX5zjdeSJSCs+4yckFgFoC60QkFegGjHK2N8Qz1pF5ns9l/JBVHzWuE5HiqprufKt8C9ipqq+5HVdeOF0t/1TVmEv0fmF4/kh3PzlOkcM+9YHHVPWeAni/PwPlVPWF8+w3Cpilqgvy+57m0rMxAuMPHhCRe/EMyKbimUVkziAiNfEMLCecKwkAqGqqM2U0PJezhjxR1Yl53HWTJYHAZWcExhgT4myMwBhjQpwlAmOMCXGWCIwxJsRZIjDGmBBnicAYY0Lc/wO679w+MN/+TwAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# Here's the data.\n", "plt.errorbar(omega, Vratio, errVratio, fmt = 'ko', markersize = 5,\\\n", " linewidth = 1.8,\\\n", " markeredgecolor = 'b',\\\n", " markerfacecolor = 'r',\\\n", " capsize = 5)\n", "plt.xlabel('Angular Frequency (rad/s)')\n", "plt.ylabel('Voltage Ratio')\n", "\n", "# Plot the best-fit line.\n", "xx = np.arange(1, 8e5, 100)\n", "plt.plot(xx, LRCFunc(xx, a_fit[0], a_fit[1], a_fit[2]), 'k-')\n", "plt.axis((0, 4.6e5, 0, 0.75));" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.8.8" } }, "nbformat": 4, "nbformat_minor": 5 }